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Abstract 

Under tidal forcing, icy satellites with subsurface oceans deform as if the surface 
were a membrane stretched around a fluid layer. ‘Membrane worlds’ is thus a fitting 
name for these bodies and membrane theory provides the perfect toolbox to predict 
tidal effects. I describe here a new membrane approach to tidal perturbations based on 
the general theory of viscoelastic-gravitational deformations of spherically symmetric 
bodies. The massive membrane approach leads to explicit formulas for viscoelastic tidal 
Love numbers which are exact in the limit of zero crust thickness. Formulas for load 
Love numbers come as a bonus. The accuracy on k2 and h2 is better than one percent if 
the crust thickness is less than five percents of the surface radius, which is probably the 
case for Europa and Titan. The new approach allows for density differences between 
crust and ocean and correctly includes crust compressibility. This last feature makes it 
more accurate than the incompressible propagator matrix method. Membrane formulas 
factorize shallow and deep interior contributions, the latter affecting Love numbers 
mainly through density stratihcation. I show that a screening effect explains why 
ocean stratification typically increases Love numbers instead of reducing them. For 
Titan, a thin and dense liquid layer at the bottom of a light ocean can raise k2 by 
more than ten percents. The membrane approach can also deal with dynamical tides 
in a non-rotating body. I show that a dynamical resonance significantly decreases the 
tilt factor and may thus lead to underestimating Europa’s crust thickness. Finally, the 
dynamical resonance increases tidal deformations and tidal heating in the crust if the 
ocean thickness is of the order of a few hundred meters. 
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1 Introduction 


Tidal Love numbers are three numbers quantifying the response of a spherically symmetric 
body to tides or to changes in rotation or orientation. Their computation is required 
for all applications in which global deformations intervene: tidal or despinning tectonics, 
tidal heating, true polar wander, tidal currents in Titan’s seas (applications discussed in 
Beuthe 2015 except for the last one, see Tokano et al. 2014]). Conversely, measuring 


Love numbers helps to constrain interior models. 

Membrane worlds refer to planetary bodies with a thin shell floating on a liquid layer 
Beuthe, 2015 . ‘Thin shell’ means here that deformations can be predicted with simple 


membrane equations instead of the more complicated thick shell theory. In practice, mem¬ 
brane theory applies to shells having a thickness less than five to ten percents of the surface 
radius. The term is thus perfectly suited to the large Galilean and Saturnian icy satellites 
for which electric, magnetic (including auroral), and gravity data point to the existence of 
a global ocean close to the surface (Table [I). Though observations are still lacking, Triton 

many 


Nimmo and Speneer 2015; Hand. 2015 


and Ceres are candidate membrane worlds 
smaller bodies could also enclose an ocean but are unlikely to satisfy the membrane as¬ 
sumption \Hussmann et al. . 2006 . In this paper, I choose Europa and Titan as case studies 
because of the available data, their potential for future missions, and their differences in 
internal structure and orbital period (Table [^. 


Table 1: Large icy satellites: some constraints on their crust thickness d (absolute and relative to 
the surface radius R) from gravity (G), magnetic (M), auroral (A), and electric (E) data. 


d (km) d/R (%) Data Reference 


Europa 

< 170 

< 11 

G 

Anderson et al. 

199 

8 


< 200 

< 13 

M 

Zimmer et al. 

2000 



< 15 

< 1 

M 

Hand and 

Chyba 

2007 

Ganymede 

150 - 330 

6-13 

A 

Saur et al. 

2015 



Gallisto 

< 300 

< 12 

M 

Zimmer et 

al. 

2 

000 


Titan 

55-80 

2-3 

E 

Beghin et al. 

2012 



In a previous paper, I obtained analytical formulas for tidal Love numbers using thin 
shell theory in the membrane limit \Beuthe 2015 . Although the method was by and 


large successful (especially regarding depth-dependent crustal rheology), it was lacking in 
some respects. First, it required that the floating shell be of the same density as the 
underlying ocean. In the membrane limit (shell of vanishing thickness), this is equivalent 
to assuming that the membrane is massless. I will thus call this method the massless 
membrane approach. Second, accurate benchmarking of the tilt factor formula revealed a 
mismatch associated with shell compressibility. Apparently, the classical equations of thin 
shell theory are not completely satisfactory regarding their dependence on compressibility. 
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Table 2: Bulk and orbital parameters of Europa and Titan. 


Parameter 

Symbol 

Europa 

Titan 

Unit 

Spin rate“ 

LO 

2.048 

0.456 

lO-’^s-i 

Surface radius*” 

R 

1560.8 

2574.76 

km 


GM 

3202.74 

8978.14 

km*^ s“^ 

Moment of inertia'”* 

Mol 

0.346 

0.341 

- 

Bulk density® 

Pb 

3013 

1881.5 

kg m“^ 

Surface gravity® 

9 

1.315 

1.354 

_9 

ms 

Dynamical parameter-! 

Quj 

4.98 

0.395 

10-“^ 


TPT, gatpllitp 


Nimmn pf 


/pfig Pi nT 


Anderson 


nj 


phjemerides (http: 
for Europa, 


2007 


2IIIII 


et al. 


Tnr Titan 


/ /Qgn jpi rta 


Mitri et at. 

2014 


1998; 


for Europa, 


computed trom i 
f computed from Eq. (20|l. 


less et al. 

2012 

.874 X i0“ 

li 


/»• 

for Titan. 


for Titan. 


kg 


-io-2\ 


For these two reasons, I develop in this paper an alternative membrane formalism, called the 
massive membrane approach, which is based on the viscoelastic-gravitational equations used 
to predict tidal deformations in thick shell theory. These equations have been extensively 
validated through their accurate prediction of the frequency spectrum of Earth normal 


modes Dahlen and Tromp . 1999 


Though technically complex, the massive membrane approach is based on two simple 
ideas. The hrst idea consists in using the viscoelastic-gravitational equations in order to 
propagate the three unknown Love numbers from the surface to the crust-ocean boundary, 
where they must satisfy two conditions called free-slip and fluid constraint. This procedure 
results in two relations between tidal Love numbers, the In — hn and kn — hn relations, which 
depend on the effective viscoelastic parameters of the crust. The second idea consists in 
factorizing the shallow interior from the deep interior: in the static limit of equilibrium 
tides, the Love numbers of the body with its viscoelastic crust are expressed in terms of 
the Love numbers of a simpler model (or fluid-crust model) in which the crust is fluid-like. 
Combining these two ideas leads to explicit formulas for Love numbers in terms of crustal 
parameters and of the deep interior structure. If tides are dynamical, fluid-crust models 
must be given up but it remains possible to derive membrane formulas for Love numbers in 
a model with an inhnitely rigid mantle. I will show that a dynamical resonance increases 
surface deformations and tidal heating in the crust as the ocean becomes shallower. 

The massive membrane approach is more than ‘yet another method’ for computing Love 
numbers. It has the interesting feature that there is an overlap, but no coincidence, between 
the domains of validity of the membrane approach and of the standard methods (Fig.[^. To 
be more clear, three successive assumptions are common when computing Love numbers. 
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Figure 1: Domain of validity of the membrane approach in comparison with other methods 
(‘homog. incomp.’ = homogeneous and incompressible). 


First, the interior structure of the undeformed body is assumed to be spherically symmetric. 
Without this assumption, Love numbers do not really make sense although the Love number 
concept is sometimes extended to flattened bodies in rotation \Wahr 


integration methods must be used if no other assumption is made [e.g. Tobie et al. , 2005 


Numerical 


Second, the static limit is often applied because numerical codes tend to diverge at tidal 
periods if the body contains a fluid layer (tides are particularly slow on Titan). For example, 
Wahr et al. 2006 , Rappaport et al. [2008] , and \Wahr et al. [2009] use a code assuming 


the static limit in all layers whereas Mitri et al. 
to the ocean. 


2014 only apply the static assumption 

Third, 


Numerical integration remains necessary in the static limit. Third, the 
interior structure is often discretized as an onion-like superposition of incompressible and 
homogeneous layers. These rather strong assumptions lead to the incompressible propagator 


matrix method [e.g. Sabadini and Vermeersen , 2004 which provides analytic solutions for 


two- or three-layer models while models with more layers are easily solved numerically (the 
propagator matrix method also exists in a dynamical and compressible version which is 
seldom used for reasons explained in Appendix [F]). The matrix method is stable when solid 
layers become fluid-like, contrary to most numerical codes. These qualities make it popular 


in planetology [e.g. Moore and Schubert\ 

2000| \Hussmann et al. 

2002 

Roberts and Nimmo 

2008 

Jara-Orue and Vermeersen ^ 

2011 

. By contrast, the membrane approach is based 


on the thin shell approximation, but it does not require incompressible and homogeneous 
layers. Compared to the propagator matrix method, the massive membrane approach 
is simultaneously more restrictive (requiring a thin shell) and more general (allowing for 
compressibility). As dynamical effects can be included in some cases, one could say the 
same with respect to codes computing static Love numbers by numerical integration. 
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2 Viscoelastic-gravitational theory 


This section reviews the basics of viscoelastic-gravitational theory that are used in the 
membrane approach. 


2.1 




functions and Love numbers 


Viscoelastic-gravitational theory describes the deformations of a self-gravitating body with 
a spherically symmetric internal structure. The deformations can result from tides, rota¬ 
tional flattening, surface loading or free oscillations due to an (earth)quake. In the standard 
formalism oi \Alterman et oL [1959 , the equations of motion and Poisson’s equation form a 
set of six differential equations of first order, the solutions of which are six radial functions 
Hi {i = 1,...,6). Following the conventions of Takeuchi and Saito [1972 , yi and 7/3 are 


( 1 ) 


scalars associated with radial and tangential displacements, respectively, 

K, «„.«*) = (ft. 93 
while 1/2 and are scalars associated with the stresses having a radial component: 

i^rr,CTre,ar^)= (^y2,yA-^,^^^U. ( 2 ) 

For tidal deformations, U is the tidal potential component of harmonic degree n at the 
surface and (r, 9, 4>) are the usual spherical coordinates (radius, colatitude, longitude). Tan¬ 
gential stresses can be expressed in terms of {yi,y^) \Takeuchi and Saito , 1972, Eq. (79)]. 
The total gravity potential perturbation inside the body is expressed as 


^ = y5U, (3) 

while 7/6 is related to the derivative of 7/5 (see Eq. (0 below). The yi are continuous 
within the body with one exception: 7/3 is discontinuous at fluid/solid interfaces because 
solid layers freely slip on fluid layers. The variables are nondimensionalized with the ocean 
density p, the surface gravity g, and the surface radius R: 

{gyi, 92/9 , gys , y^/p , 2 / 5 , Ry^i) ■ (4) 


Though there are six independent solutions, only three are regular at the center, leaving 
only three degrees of freedom. Eor tidal deformations of degree n, these regular solutions 
must be combined in order to satisfy three boundary conditions at the surface \Saito , 1974| : 


(2/2(7?) ,7/4(7?) , 7 / 6 ( 7 ?)) 



2n-h 1\ 

7 ? ;■ 


(5) 






















The boundary conditions on y 2 and 1/4 result from imposing that the surface stress has no 
radial component, whereas the condition on yg results from a boundary condition on the 
gradient of the gravity potential. For surface loading of degree n, the boundary conditions 
\Saito 


are 


1974 


,yi{R) ,yi{R)) = (^-5^ pj . 0. 


( 6 ) 


where pi, is the bulk density and the superscript L denotes surface loading. 

The response of the satellite to a tidal perturbation of degree n is parameterized by the 
radial, tangential, and gravitational tidal Love numbers: 


{hn,ln,kn) = {gyi{R), gy 2 ,{R), y^{R) - 1). (7) 

By analogy, the response of the body to surface loading of degree n is characterized by 
three load Love numbers: 


{hnJn,K) = [9v{{R) :gyhR) - 1 ) • 


( 8 ) 


In pre-spatial geodesy, measurements at the surface yielded combinations of Love num¬ 
bers \Melchior 1978; Hussmann et al. 2011| , three of which are the tilt factor 7 „ (also 
called diminishing factor), the gravimetric factor Sn, and the strain factor Sn- 


7n 

— It k^i , 

(9) 

dn 

= (n +2hn - (n +1) kn)/n, 

(10) 


= 2h„ - n (n-h 1) In ■ 

(11) 


Similar combinations 5^, can be defined for load Love numbers. Nowadays, Doppler 
and range tracking of an orbiting spacecraft and satellite laser altimetry lead to indepen¬ 
dent estimates of k 2 and /i 2 , respectively (measurements of I 2 with this technique are 
unsatisfactory). 

In this paper, I only use 6n and Sn as compact notations. By contrast, the tilt factor 
is amenable to direct physical interpretation even if surface measurements are unavailable. 


First, 72 is required to predict the height of equilibrium tides in the seas of Titan Tokano 


et al. , 2014| : it represents the reduction that must be applied to the surface tide because 


the bottom of the sea also deforms (in which case it is more appropriately called the 
diminishing factor). Second, 72 is of great interest for estimating the crust thickness of 
icy satellites with a subsurface ocean. On the one hand, 72 is known once the tidal Love 
numbers k 2 and /i 2 have been determined from gravity and altimetry data. On the other, 
72 is approximately proportional to the crust thickness dii d/R 1 \Wahr et al 


2006 


Measuring the tilt factor is thus one way of estimating the crust thickness. Compared to 
the Love numbers taken separately, the tilt factor has the advantage that it is less affected 
by uncertainties about the interior density structure. 
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2.2 Viscoelastic-gravitational equations 


Consider a body whose properties vary only with the radius: density Prir), unperturbed 
gravitational acceleration Prir), and Lame constants A(r) and p{r) (the latter is also called 
shear modulus). In order to facilitate comparisons with thin shell theory, I use Poisson’s 
ratio u{r) instead of the first Lame constant A(r) (see Table [^. For later use, I define the 
compressibility factor x as: 


X = 


1 - 

1 - i' ■ 


( 12 ) 


For common materials, the compressibility factor ranges from y = 0 (if = 1/2 i.e. 
incompressible limit) to y = 1 (if = 0). For elastic ice, z/ ~ 1/3 so that y ~ 1/2. Table 
summarizes the relations between the various compressibility parameters (A, K, v, x) and 
common pairs of elastic constants. 


Table 3: Various compressibility parameters 



(h, A) 

iP,K) 

{P. > 2 ) 

A 

A 

k-Ip 

2/iz/ 

l-2i^ 

K 

A + 

K 

3 l-2v^ 

y 

A 

3K-2IJ. 

y 

2(A-|-^) 

QK+2n 

X 

2m 

6/j. 

l-2u 

A-1-2^ 

3if+4M 

1-12 


For deformations of harmonic degree n and angular frequency cj, the viscoelastic- 
gravitational equations (Eq. (82) of Takeuchi and Saito [1972]) can be written as 


y'l 


y2 


y's 


y4 


2/5 


2/2 - ^ (1 - X) ( 22/1 - n (n + 1) 1 / 3 ), 
2p r 


2 2 1 / 1 + u 

= -UJ PrVl - - xy2 + ^ \ 


- prPrT (2yi - n (n + 1 ) 1 / 3 ) 


n(n + 1) 

H- VA- Pr\ ye 

r 

= -2/4 + ^ (2/3 - 2/1), 

p r 

2 1 / N 2 

= -w /ir 2/3 - - (1 - X) 2/2 - ^ At 

- - 2/4 - — (2/5 - 9 ryi ) , 
r r 

. „ 7i + 1 

= y & + 47 rGpr- 2/1- 2/5 , 


1 - V 

n+1 2 

- 2/5 + - 2/1 

r r 


1 + 12 
1-12 


2/1 


Xn + 1 + T2 
1 - 12 


2/3 


2/6 = 


n - 1 AirGpr . , , w x 

-2/6 H- [n + 1) [yi - nys ), 


(13) 

(14) 

(15) 


r 


r 
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(16) 

(17) 

(18) 



























where the prime denotes a derivative with respect to radius. The radial dependence of the 
functions yi and of the parameters {pr, gr, is implicit. In Eq. (16), Xn is defined by 


Xn = {n- 1) (n + 2), (19) 

which is the degree-n eigenvalue of the operator —(A+ 2), A being the spherical Laplacian. 


2.3 Static versus dynamic 


Dynamical (or inertial) terms are the terms of Eqs. (14) and (16) that are proportional to 
Lo"^ and which result from the acceleration term in the equations of motion. Compared to 
seismic perturbations, tidal deformations are slow so that dynamical terms are expected 
to be small with respect to other terms. This assertion can be verified by nondimen- 
sionalizing Eq. (14) with Eq. 0 and comparing the different terms which are of three 


types: dynamical, viscoelastic, and gravitational. Dynamical terms are parameterized by 
the dimensionless quantity q^j'. 

quj = 




9 


( 20 ) 


Eccentricity and obliquity tides on synchronously rotating satellites have a frequency oj 
equal to the spin rate of the body. In that case, is equal to the ratio of the centrifugal 
acceleration at the surface to the gravitational acceleration, which is a parameter (denoted 
q or m) commonly used in the theory of equilibrium figures [e.g. Schubert et al. 2004| . 


Viscoelastic terms are parameterized by the nondimensional shear modulus, 


f-^ice — 


l^ice 

PiceQ^ 


( 21 ) 


where pice 3.5 GPa and pice ~ 1000 kg/m^ while ( 5 , R) are given in Table Finally, the 
nondimensional quantity associated with gravity terms is 1 . 

In the crust, dynamical terms are small with respect to viscoelastic terms if <C Pice 
whereas they are small with respect to gravitational terms if ^ 1 • The former constraint 
means that the tidal velocity ojR is small with respect to the seismic S-wave velocity 
Pice!Pice- The latter means that the tidal frequency u is small with respect to the free 
oscillation frequency of an incompressible homogeneous fluid sphere of radius R and surface 
gravity g (Eq. (124) with z = 0 and ^ = 1; or Eq. (8.188) of Dahlen and Tromp [1999] ). 
For large satellites, the two constraints are similar because pice ~ 1.7 and 1.0 for Europa 
and Titan, respectively (for smaller bodies, the viscoelastic constraint is weaker). Since 
go; 1 for Europa and Titan (Table [^, it is an excellent approximation to neglect these 
terms in the crust as long as there is a lithosphere, i.e. the crust is mechanically rigid close 
to the surface. The same argument holds for the mantle and the core (if it is solid) by 
substituting pice and pice with the appropriate shear modulus and density of the layer. 
Dynamical effects, however, can be significant in a liquid layer (see Section]^. 
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In a given layer, the static limit consists in setting 


oj = 0 


(22) 


in the viscoelastic-gravitational equations. The static limit is typically applied to the whole 
body 


[e.g. 


Wahr et al.. 2006 or only to the ocean in order to stabilize the solution Mitri 


et al., 2014 . In Section sjT will do the opposite by taking the static limit in the solid layers 


but not in the ocean. This choice has the advantage of preserving the dominant dynamical 
effect while considerably simplifying the solution. 


2.4 Fluid layer 

The case of a fluid layer (here the subsurface ocean) deserves particular attention for two 
reasons. First, free slip between ffuid and solid layers leads to new boundary conditions 
at their interfaces. In particular, free slip at the crust-ocean boundary (of radius i?^) is a 
crucial condition in the membrane approach; 

y4{Re) = 0. (23) 

Second, displacements within the fluid are undetermined in the static limit, in which case 
only gravity variables can be propagated through the ffuid. A fluid layer is characterized 
by vanishing shear stress and shear modulus (?/4 = 2/4 = 0 and /x = 0 ) which implies 
that u = 1/2 (from its definition in Table |^. Within the fluid, the fourth viscoelastic- 
gravitational equation (Eq. ([I6|)) becomes 

y 2 = p {gryi - 2 / 5 ) - 2 / 3 , (24) 


where p is the ffuid density. 

In the static limit [oj = 0), displacement and stress become indeterminate inside the 
ffuid layer (though the radial variables 2/1 and 2/2 are well-defined at ffuid/solid interfaces). 
Nevertheless 2/1 and 2/2 are constrained by Eq. (24) with a; = 0: 


2/2 = /5 {gryi - 2/5) • 


(25) 


In spite of this indeterminacy, knowing the potential and its derivative in the ffuid is 
sufficient to solve the viscoelastic-gravitational problem. As 2/6 depends on the displacement 
2 / 1 , it is replaced by the variable 2/7 which depends only on the gravitational potential and 
is everywhere continuous: 


47rG 

2/7 = 2/6 4-2/2 • 

gr 


(26) 


The equation for y^ is replaced by an equation for 2/7 involving only the variables 2/5 and 
2/7 \Saito I974| . Within the ffuid layer, the gravity potential (in the static limit) is thus a 
superposition of two general solutions decoupled from stress and strain. 
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If the fluid layer reaches the surface (surface ocean or quasi-fluid crust), the boundary 
condition 7/2 (-R) = 0 (Eq. Q) applied to the static fluid constraint (Eq. ([^) leads to a 
relation between Love numbers: 

K + l = K, (27) 

where the superscript ° indicates that the surface layer of the body behaves as a fluid. 
This equation means that tides are in hydrostatic equilibrium; the surface of the satellite 
coincides with the geoid. 


3 The crust as a membrane 

In this section, I apply the membrane approximation to the viscoelastic-gravitational equa¬ 
tions. As a result, crustal rheology will be described by three (if tides) or four (if surface 
loading) effective viscoelastic parameters. If there is no bulk dissipation, only two of them 
are independent. 


3.1 Principle 


As mentioned in the introduction, the basic idea of the membrane approach consists in 
propagating the yi variables from the surface, where they are either fixed by the boundary 
conditions or parameterized in terms of Love numbers, to the crust-ocean boundary where 


two constraints hold: the free-slip condition and the fluid constraint (Eqs. (23)-(24)). The 


idea is depicted in Fig.[^ If the crust thickness d is small with respect to the surface radius 
R, all equations can be expanded in the small parameter e: 


d 

^~R' 


The radius of the bottom of the crust (or crust-ocean boundary) is denoted 

Re = Ril - e). 


(28) 


(29) 


I will evaluate the variables yi at the bottom of the crust in the thin crust limit, that is at 
first order in £ (denoted 0(e)). 

The value of y* at the bottom of the crust is related to Its surface value by 


fR 

yi(Rs) = yi(R) - / y'i(r)dr. 

J Re 


(30) 


The derivatives y( are given by Eqs. (13)-(18). The integral over the crust thickness being 


proportional to e, I only need to keep terms In y( that are 0(1) (i.e. zeroth order in e). 


Thus I can make the following approximations in the right-hand sides of Eqs. (13)-(18): 


1. Evaluate all quantities at r = R, except the density and the viscoelastic parameters. 
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hnl kn'l 1,1 

ym 



+ Fluid constraint 
^ In—hn and relations 

Figure 2: First basic idea of the membrane approach: propagate the yi variables from the surface 
(radius R), where they are either fixed by the boundary conditions or parameterized in terms of 
Love numbers, to the crust-ocean boundary (radius R — d) where the free-slip condition and the 
fluid constraint must be satisfied. The results are two relations between Love numbers. 


2. Apply the surface boundary conditions for tides or surface loading (Eq. Q or ([^). 
Express yi{R), VsiR), and y^{R) in terms of Love numbers (Eq. 0 or ([^). Com¬ 
binations of Love numbers can be compactly written with the factors 7 n, (5n, and Sn 
(Eqs. (§-(0) or with the corresponding surface loading factors ( 7 ^^, 5'^^, s',). 


3. Static limit in the crust: set a; 


0 in Eqs. (14) and (16). 


3.2 Density and gravity 


In the membrane approximation (Eq. (30)), density terms are integrated over the crust 
thickness. What remains is the mean density of the crust: 



(31) 


In general, the ocean density increases with depth. The density of the top layer of the 
ocean, in contact with the crust, is denoted p. The density contrast at the crust-ocean 
boundary is denoted 

Sp = p-p. (32) 

The bulk density (or mean density of the whole body) is denoted pb- The ocean-to-bulk 
and crust-to-bulk density ratios are defined by 




P_ ^ 
Pb ’ PbJ 


(33) 
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The surface gravity (g) and the gravity at the crust-ocean boundary (g^) are given by 


dvr 

g = Y ’ 


g, = {l + e{2-3i))g. 


(34) 

(35) 


The second equation is valid at 0(e). It results from g^ = {A-k/^)G peRe, where pe is the 


mean density of the body without its crust. At 0(e), Pe = Pb + 3{ph — p)e so that Eq. (35) 
follows. 


3.3 Effective viscoelastic parameters 

As rheology depends on depth, viscoelastic parameters cannot be considered as constant 


when integrating y[ over the crust thickness (Eq. (30)). In Eqs. (13)-(18), the parameters 


p, V and X appear in y\ in various combinations. For tidal deformations, two of them (x!P 
and l//r) are eliminated by the zero boundary conditions on y^ and 1 / 4 . It is thus only 
necessary to integrate on x and on the parameters (p, q) defined by 


(p, g) = 


P 


The inverse relations read 


pv 


1 - 1 - 


(p,v)= [p-q,- 

\ p 


(36) 


(37) 


Integrating x and (p, q) on depth yields three effective parameters (denoted with a bar): 


X = 


pu 


dr. 


(38) 


The effective shear modulus p and the effective Poisson’s ratio v are then defined by 


relations similar to Eq. (37): 


{p,v)= [p-qX 

p 


which can be inverted as 


(p, q) = 


P 


pv 


1-p’ l-v 


(39) 


(40) 


In other words, I have defined the effective parameters so that integrating over the crust 
thickness amounts to put a ‘bar’ on the parameters: 


X 


P 


pu 

1 — U 1 — ly 




(41) 
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In terms of (/U, u), the effective parameters u) are given by 




v = 




1 - V 


dr 


-1 


1 - V 


dr. 


(42) 

(43) 


The interpretation of x and fi is straightforward; x is the mean of y and p, is the mean of 
/i. The definition of is a bit more complicated: i> is the weighted mean of i' with weight 
p. Remarkably, p and i> are identical to the effective shear modulus and effective Poisson’s 
ratio of the massless membrane approach \Beuthe 2015 . The effective compressibility 


factor T) however, does not appear in the massless membrane approach. 

In general, the relation between x and (p, q) is not linear so that x cannot be written as 
a combination of p and q. The assumption of zero bulk dissipation, however, imposes the 
condition that K = Ke (the subscript E stands for ‘elastic’). Since K = {2/3){p + q)/x 
(see Table , the condition K = Ke is equivalent to 

X ^ p + q 
Xe Pe + Qe 


(44) 


If the elastic parameters {pe,ee) are constant with depth, integrating this equation on 
depth yields 


X = 


XE 


(p + q) 


PE + QE 

1 — 2h'E P + p 
1 + I^E PE - P 


(45) 


This relation shows two things: 

• there are only two independent effective parameters: p and P. 

• X varies with viscosity approximately like p/2pE (if ee ~ 1/3) because P is much 
less sensitive to viscosity than p (Fig. [^. 

If bulk dissipation is not zero or if {pE, ee) vary with depth, x is an independent parameter 
which must be computed with Eq. 


(38) 


In membrane equations, the effective shear modulus is always multiplied by the mem¬ 
brane thickness d. It is thus often replaced by the extensional rigidity = 2pd which 
relates integrated plane stresses to plane strains (Eq. (5) of Beuthe 2015| ). The effective 
shear modulus and the extensional rigidity are nondimensionalized as follows: 


P = 




P 

pgR ’ 

2p 

1-p' 


(46) 

(47) 
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Computing load Love numbers requires integrating on x//U (Eq. p^) or, equivalently, 
on the nondimensional parameter = x//i: 

Xm = 


The second line holds if there is no bulk dissipation (same argument as above) and shows 
that Xfi is not an independent parameter under that assumption. 


1 - ‘2.1'e pgR 
1 + ue pe 


(3 - 2x) 


(48) 


3.4 Example: stagnant lid regime 


As an illustration, suppose that the crust is made of two uniform layers: the top layer is elas¬ 
tic (no viscous effects) while the bottom layer has a rheology of Maxwell type (Appendix [A]). 
If dtop is the thickness of the top layer, its relative thickness is denoted / = dtop/d. This 
two-layer structure can represent the rheology of either a conductive crust \Ojakangas and 
Stevenson. 1989| o r a conductive/convective crust in the stagnant lid regime [e.g. 


lUSS- 


mann et al.. 2002|. In the former case, the crust is nearly elastic throughout and the two 


layers are not distinguished. In the latter case, convection occurs under a conductive lid: 
the top layer is conductive and elastic whereas the bottom layer is convective and vis¬ 
coelastic. Fig. shows the effective viscoelastic parameters in terms of the dimensionless 
number 5 which is inversely proportional to the viscosity of the bottom layer (Eq. (A.3)). 
On the x-axis, values from left to right correspond to a bottom layer which is elastic-like 
(6 < 0.1), critical (0.1 < 5 < 10), and fluid-like (6 > 10). If the bottom layer is fluid-like, 
the absolute value of the effective shear modulus drops to a fraction / of its elastic value: 
\JA ~ fpE or, equivalently, \p\d ~ PEdtop- In that case, the crust response is determined 
by the elastic top layer which acts as a lithosphere. More generally, it makes sense to define 
the lithospheric thickness by 

diitho = — d. 

PE 


(49) 


The classification elastic/critical/ffuid-like is based on the dependence of p on viscosity 
and would be shifted to higher values of 6 were it based on the behaviour of P (compare 
panels A and C in Fig. [^. If the top layer is purely elastic, the critical transitions occur at 


• 5 = 1: middle point of Re{p) and maximum of Im{p), 

• 5 = 5': middle point of Re{x) and maximum of Im{x)-, 

• 5 = 5": maximum of Re{p) and zero-crossing of Im(p), 


where 5' = 3(1 — z^e)/(1 -|- ue) and 5" = \/5'jf. 
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Figure 3: Effective viscoelastic parameters as functions of inverse viscosity for a conduc¬ 
tive/convective icy crust: (A) shear modulus (normalized by its elastic value), (B) compressibility 
factor, (C) Poisson’s ratio. The crust is divided into a top elastic layer and a bottom viscoelastic 
layer with / being the relative thickness of the top layer. The elastic Poisson’s ratio is = 0.33. 
Black (resp. gray) curves correspond to a top layer that makes 40% (resp. 10%) of the total crust 
thickness. Solid (resp. dashed) curves show the real part (resp. imaginary part). Vertical dotted 
lines indicate the positions of the critical transitions. See Section 


3.3 


for details. 


3.5 Propagating yi functions through the crust 


The yi functions can now be propagated from the top to the bottom of the crust. Using 
Eqs. (38)-(40) and the approximations of Section 3.1, I integrate Eqs. (13)-(18) and insert 
them into Eq. (30). Eor tidal deformations, I get the following nondimensional equations: 


gyiiRe) 

= hn + e{l- x)sn, 


(50) 

-2/2 (-Re) 

= ~ ^(1 + p) bex ~ 

P\ P r 

^ — 1 T ^ ; 

(51) 

p 

Pj P 

gysiRe) 

— (1 Ifi -\r £ hfi , 


(52) 

^ 2/4 (-Re) 

— Dex ^ ~ 

{Xn + 1 + P) In) + e - 7n ) 

(53) 

p 

^ P 

ybRe) 

= (1 + (n + 1 ) 6) {kn 

“b 1) — 3 e ^ hn — (2?r -b 1) £ , 

(54) 

Ry&{Re) 

= (2n + l)(l-(n- 

1)£) - 3e^((n -l)hn + Sn) ■ 

(55) 


I will also need the auxiliary gravity variable y^ at the crust-ocean boundary. Substituting 
Eqs. and (55) into Eq. (26), I get 


Ry^{Re) = (2n -b 1) (1 - (n - 1) e) - 3e ^ ((n - 1) hn - n5n) - 3 ^ b^x b + i^) Sn ■ (56) 
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Finally, the quantity grUi — 2/5 is of particular interest because it appears in the fluid 
constraint (Eq. (24)). Combining Eqs. (35), (50), and ([54)), I can write 


geVliRe) - VbiRe) = -Jn + £ n6n + S {I - x) Sn ■ (57) 

4 Relations between Love numbers 

In Section I propagated the yi functions from the surface to the crust-ocean boundary 
(Eqs. (@-(1^). I now examine the constraints imposed by the special kind of coupling 
between crust and ocean: the ocean exerts a radial push on the crust (continuity of the 
radial stress) but no lateral traction (shear stress vanishes, free slip occurs). Table [^ 
summarizes the parameters relevant to the shallow interior which will be used repeatedly 
in the rest of the paper. 

Table 4: Parameters relevant to the shallow interior. 



Parameter 

Symbol 

Size 

Surface radius 

R 


Crust thickness 

d 


Relative crust thickness 

e 


Radius of crust-ocean boundary 

Re 

Density 

Mean crust density 

P 


Ocean density (top layer) 

P 


Crust-ocean density contrast 

Sp 


Bulk density 

Pb 


Ocean-to-bulk density ratio 



Crust-to-bulk density ratio 


Gravity 

Surface gravity 

9 


Gravity at crust-ocean boundary 

9e 

Elasticity 

Effective shear modulus 

P 


Effective Poisson’s ratio 

D 


Effective compressibility factor 

X 


Effective shear modulus (nondim.) 

p 


Extensional rigidity (nondim.) 

D ex 

Varia 

Eigenvalue of — (A -|- 2) 

Xn 


Dynamical parameter 

quj 


Eq. 


(29) 

(31) 


(32) 


(^) 

(^) 

(M) 

(^) 

(42) 

(43) 

(45) 

(46) 

(47) 
(W) 
( 20 ) 
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4.1 Magnitude of tilt factor 

Let us start with the radial stress coupling at the crust-ocean boundary. At the top of the 


ocean, the fluid constraint (Eq. (24)) reads 


^y2{Re) = geViiRs) - VbiRs) - ^^^-^gyi{Re 


(58) 


As the crust freely slips on the ocean, the function is discontinuous at the crust-ocean 
boundary: the superscript F indicates that must be evaluated on the fluid side of the 
boundary. I will assume for the moment that y^(Re) and 7/i(i?£) are of the same order of 
magnitude (this assumption is discussed in Section]^. Since tides are slow (Section |2.3[ ), 
the last term in Eq. (58) is of 0(e) or smaller. 

Contrary to the function y^, the functions yi, y 2 and y^ are continuous at the crust- 


ocean boundary. I can thus substitute Eqs. (51) and (57) into Eq. (58). The result is a 


relation between Love numbers, or rather their linear combinations 'jm and Sn- 


e 

7n = ( (1 + 1^) Dex - £-ex 

p 


6p 


-e^n5n-quj gy^ (Re) ■ 


(59) 


All terms in the right-hand side are of 0{e) except maybe the first one depending on the 
extensional rigidity. Two cases must be considered: 

• the crust is soft (lHexI ^ !)• Love numbers are large and the tilt factor is small. 


(hri) I'm kni ^rii Sr, 


In 


0 ( 1 ), 

0(e), 


(60) 


which also means that 


the crust is hard {\D, 
factor is close to one, 


kn F 1 — hn 0(e). (61) 

exi ^ 1 and DexSn ~ 0(1)). Love numbers are small and the tilt 


(h-ri) kn, 6n, Sr 


In 


< 1 , 

~ 1 . 


(62) 


Erom now on, I assume that the crust is soft so that Eq. (60) holds. In this way, I will 
obtain formulas for kn and that are valid at 0(e). If the crust is hard, the extensional 
rigidity Df,x must be considered as a parameter of 0(1). In that case, the expansions 
of the viscoelastic-gravitational equations (Eqs. (50)-(55)) are not complete at 0(e): one 
should indeed include terms like Dex^, which can only be obtained by formally expanding 
the viscoelastic-gravitational equations to O(e^). All is not lost, however: up to 0(1), 
the formulas for a soft crust also apply to a body with a hard crust. This means that, if 
the crust is hard, only the dominant contribution of the crustal rigidity must taken into 
account, while density terms and other small corrections must be neglected for consistency. 
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4.2 Relation between and hn 


Imposing the free-slip condition yi{Re) = 0 on Eq. (53) yields In in terms of {hn,jn) at 
the crust-ocean boundary: 


R. — 


l + u 


Xn + l + l' 


, , 1 

- hn -\ -— 


1 - i? 


p 2jl Xn + I + l' 


-In- 


(63) 


If the crust is soft, the first and second terms in the right-hand side are of 0(1) and 0(e), 
respectively, because 7 n ~ 0{s) (Eq. (60)). Going back to Eq. (53), one sees that the 
second term of Eq. (63) comes from a term of O(e^) and should in principle be neglected 


because the equations were not expanded beyond 0(e). This term, however, is not small in 
the limit of a quasi-fluid crust because of the prefactor l/p. I will exclude this possibility 
by imposing the following lithospheric condition (a more precise form will be given later): 


P > In/hr. 


(64) 


This constraint holds as long as there is a lithosphere, i.e. the upper part of the crust has 
nonzero rigidity. Under that reasonable assumption, the free-slip condition yields a relation 
between the displacement Love numbers, or In — hn relation: 


In — 


l + P 


hn 


ill — 41 ? 

Xn + l + n 


(65) 


which is valid at 0(1), i.e. in the limit of zero crust thickness. The same relation holds for 


a hard crust because the second term of Eq. (63) is smaller than the first one by a factor 
e (hn ~ 0(l/fie) and 7 „ ~ 1). 

The In — hn relation coincides with the one derived in the massless membrane approach 
(Appendix [B]). Einite thickness corrections depend not only on the thickness and density 
of the crust but also on its rheology \Beuthe 


2015, Fig. 4]. If the crust does not convect. 


it is nearly elastic and finite thickness corrections can be estimated with the homogeneous 
crust model of Appendix If tides are static, the I 2 — h 2 relation for this model, up to 
order 0(e), can be read from Table 11 


l^ = l(l-^e]h2. 
11 V 33 ' 


( 66 ) 


The factor 3/11 results from (x 2 -|- 1 -|- z^)/(5 -|- u) with n = 1/2. 
approximation if e < 0.3 \Beuthe , [2015 Fig. 13]. 


Eq. (66) is a good 


4.3 Relation between k„ and hr 


In Sections 4.1 and 4.2, I obtained two independent relations between Love numbers: the 
In — ^n — Sn and In — hn relations which are of 0(e) and 0(1), respectively. I will now 
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combine them so as to relate kn and h. 

Eq. (59) can be evaluated at 0{e) by knowing s- 


(65) into Eqs. ([^-(10), I get 


at 0{e). If the crust is soft, the right-hand side of 
and 5n at 0(1). Inserting Eqs. (61) and 


Sn — 


Xn (1 - 


Xn 


nS„. = 


+ l + iy 
= 2n + 1 — {n — 1) h 


hn + 0(e), 

-I- 0{e). 


(67) 

( 68 ) 


Note that the lithospheric condition, Eq. ( |64[ ), is required for the first approximation to 
be valid whereas the second approximation depends on the assumption of a soft crust. If 


the crust is hard, Eq. (68) has an additional term — (n -|- l)7n in the right-hand side. This 


modification does not matter because the right-hand side of Eq. (59) is then of 0(1) and 
all terms of 0(e), including the one depending on n6n, are neglected in that case. 


The substitution of Eqs. (67)-(68) into Eq. (59) yields a simple relation between the 


tilt factor and the radial Love number: 


6p 


7 n = Arhn- {2n + l) —e . 

P 


(69) 


In the right-hand side, the last term is the major density correction while the factor At in 
the first term is the sum of four contributions: 


At — A -|- -|- Ap Ai^ . 


(70) 


The first (and generally dominant) contribution is the membrane spring constant A, which 
vanishes if the crust is fluid-like: 


A = fpjle. 


( 71 ) 


The dimensionless coefficient fp is defined in Table together with the coefficients and 
fp appearing below. The name ‘membrane spring constant’ comes from the observation 
that the membrane radial response follows Hooke’s law (Appendix [B|). 

The second contribution to At is the compressibility correction A^, which vanishes if 
the crust is incompressible: 

^x = fxXe. (72) 


The third one is the minor density correction Ap, which vanishes if there is no density 
contrast at the crust-ocean boundary: 


Ap = fp^e. (73) 

P 

For tides of degree two, the qualification ‘minor’ is justified because this term is about ten 
times smaller than the major density correction (/p ~ 0.5 if n = 2). The fourth one is the 
dynamical correction A^^, which vanishes in the static limit: 


Afjj — Qijj 


vim ^ 

yi{Re) ' 


(74) 
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The ratio {Re)/yiiRe) is unknown at this stage but is supposed to be of order unity 
so that Aoj ~ (Section 4.1). Dynamical corrections are discussed in more detail in 

Section HI 


The tilt factor formula (Eq. (69)) can be written as a relation between Love numbers, 
or kr, — hr, relation: 


+ 1 — (IT A't’) hn — {‘2,n + 1) — £ . 

P 


(75) 

When the crust becomes fluid-like, Eq. ( |75[ ) does not tend in the static limit to the hy¬ 
drostatic — /i° relation (Eq. (27)) because of the terms proportional to 5p. This is not 
surprising because the derivation of the kn - hn relation requires the lithospheric condition. 


Eq. (64). This condition can be reformulated with the help of Eq. (69): 


Sp 

^ —e, 

P 


(76) 


How does the above kn — hn relation compare to the one obtained with the massless 
membrane approach \Beuthe 2015 ? In the latter approach, the crust and ocean have the 
same density {5p = 0) and tides are static (A^^ = 0). With these assumptions, Eq. (75) 
becomes fc„-|-l = (1-|-A-|- A^^,) hn- This relation differs from the massless membrane 
relation by the compressibility factor A^^, which accounts for the mismatch found by Beuthe 


2015 when benchmarking the tilt factor formula. The differences between the massive and 


massless approaches are analyzed in more detail in Appendix [B) 


Table 5: Dimensionless coefficients appearing in the kn — hn and k!^ — h'^ relations (Eqs. (71 )-(73l 
and Eq. (79l). 



arbitrary n 

n = 2 

u 

2Xn{l + v) 

8(1+P) 

Xti + 1 + P 

5+P 

fx 

Xr^{l — y) 

4(1-P) 



/p 

(n—l)(n^ —3+(n+3)L') 

1+5P 

Xn-\-l-\-i> 

S+b' 

/; 

{n^ — l)(n+l —i/) 

3(3-P) 

Xn + l+i> 

S+P 

fxp 

n(n+l)(l — p) 

3(1-P) 

2(£c„ + 1 +P) 

S+b' 


4.4 Load Love numbers 

Similarly to tidal Love numbers, relations between load Love numbers can be obtained 
by propagating the yi functions through the crust and applying the fluid constraint plus 
the free-slip condition at the crust-ocean boundary. The difference is that the nonzero 
boundary condition on y 2 (Eq. ([^) introduces new terms in the equations of propagation. 
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The I'n — h'^ relation reads 

l'r,= _ h'^ + 

Xn + 1 + 1 ^ 


2n + ll 1 — V /_ 6p 

3^ 2/i Xn + 1 + V p 


(77) 


The second term in the right-hand side results from the nonzero radial stress at the surface. 


As in Eq. (63), the term 1/p diverges if there is no lithosphere. The k'^ — h'^ relation reads 

2n + 1 


t; + l-(l+AT)/.;-(2n+l):^. + (l + V) 


3 ? 


(78) 


This relation has the same form as the kn — relation except for the term proportional 
to (2n + l)/(3^) which results from the nonvanishing radial stress at the surface. The 
dimensionless number ijj gathers new compressibility and density corrections of 0(e): 




1 


V’ = -/x X + /p — + X X/2 + /xp w X + — 




(79) 


where the coefficients /^, f^, and fy^p are defined by Table while Xp is given by Eq. (48). 

Now the gravitational load Love number is related to tidal Love numbers by the Saito- 
Molodensky relation [Molodensky 1977; Saito, 1978, Lambeck, 1980 


k' — h — h 
= 7n - 1 • 


(80) 


Combining this equation with the tilt factor formula (Eq. (69)) and the k'^ — h'^ relation 


(Eq. (78)), I can express the radial load Love number in terms of the radial tidal Love 
number: 

At/ in - (1 + V’) ^ 

1 + At-’ 


h' = 


where ifj is defined by Eq. (79). The computation of load Love numbers is thus reduced to 


the evaluation of tidal Love numbers. 


4.5 Rigid mantle model 

In this section, I check the validity of the kn — hn relation against a semi-analytical model 
derived with thick shell theory. Consider static tides of degree two acting on a body 
described by the rigid mantle model: the mantle is infinitely rigid and the ocean is homo¬ 
geneous and incompressible. Wahr et al. [2006 solved the degree-two elastic-gravitational 
equations of this model assuming that the crust is incompressible and of uniform elasticity. 
Since the resulting analytical formulas are very complicated, they expand them at 0{e) so 
as to obtain simple formulas for Love numbers. Finally, they fit corrections due to crust 
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compressibility with a numerical code. I will now check that membrane formulas reproduce 
these results if the crust is homogeneous, i.e. in the limit y) = (/r, x). 

I start by evaluating the tilt factor for an incompressible crust {u = 1/2) before es¬ 
timating compressibility corrections. At 0(1), /12 is given by the formula for a two-layer 
body with rigid mantle and surface ocean: = 5/(5 — 3^) (Eq. (C.l) in which —)• ^). 


Substituting this formula in the right-hand side of Eq. (69) and setting = 1/2 yields 


inc 

72 


_ /24 




Vll pgR \ll p 


f^2r ^ ) 


(82) 


where the superscript ‘inc’ denotes that the crust is incompressible. This formula coincides 
with Eq. (7) of Wahr et al. [2006 


Next, I express corrections due to crust compressibility in terms of the parameter p/K 
which is related to Poisson’s ratio v (see Table [^. In particular, v = 1/3 is equivalent to 
p/K = 3/8. Using the expansion v ~ 1/2 — p/2K, I get 


__ inc I ^ 

72 — 72 ^ ^ ' 


where 




8 + 11 + 6 ^ . 
pgR p ' 


(83) 


(84) 


The terms within brackets result from A, Aj^, and Ap, respectively. Eor Europa, Wahr 


et al. [2006] assume that ^ 1/3, p/{pgR) r\j 1 and 5p/p = —0.08. These parameters 

yield 72 "''’ ~ 3.1e and C ~ —1.5, i.e. a correction of 18%. This result is close to the 
correction of C = —1.4 that Wahr et al. [2006 found by fitting Eq. (83) to their numerical 
model (see their Eq. (12)). Eor Titan and Ganymede, the correction can be smaller or 
larger, depending on the choice of p and 5p, while it goes down to 9% for small bodies 
with p/{pgR) ^ 1. Corrections due to crust compressibility are larger than the error 
due to the membrane approximation, unless the crust is very thick (Section]^. They are 
also larger than corrections due to core-mantle elasticity (Section]^. If the crust is thin, 
the membrane approach is thus more accurate than the incompressible propagator matrix 
method. 


5 Static Love numbers 

5.1 Outline of the method 

The In — hn and kn — hn relations of Section depend explicitly on the properties of the 
shallow interior (crust and top layer of the ocean). By contrast, the deep interior structure 
appears only implicitly through the Love numbers themselves. This separation between 
the contributions of the shallow and deep interior can be taken further with the second 
basic idea of the membrane approach (Eig. [^: determine the Love numbers of the body 
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Physical model 



Fluid-crust model Massive membrane 



Figure 4: Second basic idea of the membrane approach (for static tides): decomposition of the 
physical model into a fluid-crust model to which is added the contribution of the crust (‘addition’ 
should be taken in a figurative sense). Darker (resp. lighter) shades represent solid (resp. liquid) 
layers. See Section 15.21 for details. 


with its viscoelastic crust (or physical model) in terms of the Love numbers of a simpler 
model (or fluid-crust model) in which the crust behaves as a fluid \Beuthe 2015 . In a 
sense, this method factorizes the effect of the membrane from the deep interior. I will here 
extend the method of Beuthe [2015 to a thin crust of finite thickness with a density jump 
at the crust-ocean boundary. 

The static limit is a crucial requirement of the factorization between membrane and 
deep interior. In this limit, gravity decouples from stress and strain within the fluid layer, 
so that one can solve for the gravity variables {y^-,yr) independently of the displacement 
and stress variables (yi,..., yfl) \Saito 1974 . The method consists in solving a linear system 
of two equations consisting of 

• a constraint on the gravity variables at the crust-ocean boundary, 

• a scaling relation between the gravity solutions of the physical and fluid-crust models. 


Once kn is determined, radial and tangential Love numbers result from the kn — h^ and 
In — hn relations. 


5.2 Fluid-crust model 

The fluid-crust model is defined as having the same internal structure as the body to be 
modeled (called the physical model), except that the crust is fluid and of density p°. There 
are two obvious choices for the fluid-crust density p°: it is equal either to the original crust 
density {p° = p), oi to the density of the top layer of the ocean {p° = p). The latter choice 
makes it simpler to compute the Love numbers of the fluid-crust model, but it changes the 
bulk density of the body from pi, to p^. Conservation of mass yields at 0(e) 

Pb = pl + 3 5p° e , (85) 


26 

















where 


dp° = p-p° 


( 86 ) 


In analogy with Eq. (|33|), the ocean-to-bulk and crust-to-bulk density ratios for the fluid- 

(87) 




crust model are denoted 


Following the notation of Section \2A\ I define 

• Hi as the solutions of the physical model: viscoelastic crust (A / 0) of density p. 

• y° as the solutions of the fluid-crust model: fluid crust (A = 0) of density p°. 

I adopt the same conventions for the corresponding Love numbers. For example, 

kn + l = y^{R) , 

+ 1 = 2/5(7^) ■ 

Simple fluid-crust models are given in Appendix 


( 88 ) 


5.3 Gravity at the crust-ocean boundary 

First, I relate the gravity perturbation at the crust-ocean boundary to its value at the 
surface (equal to kn)'- Substituting Eq. ([M|) into Eq. (54), I get 


UbiRe) = (1 -h (n -b 1 - 3 ^) e) (fen -h 1) - (2n -h 1) e . (89) 

In the fluid-crust limit, this equation reads 

2/5 (T^e) = (l-|-(u-|-l — 3 e) (fe° -b 1) — (2n -b 1) e . (90) 

Next, I relate the auxiliary gravity variable y-j at the crust-ocean boundary to 2/5 at the 
surface. Starting with Eq. (56), I express the right-hand side of this equation in terms of 
hn with the help of Eqs. (67)-(|6^. Then I eliminate in favor of kn with the kn — hn 
relation (Eq. Q). Terms of O(e^) are neglected except in the term in A/(l -b A) so that 
the limit of large A is well-behaved. At the crust-ocean boundary, the gravity variable yj 
is thus given by 


RyiiRe) + 3C 


-b 2 (n - 1) ) (fen -b 1) = (2n-b 1) (1 - e (n - 1 - 3^)) . (91) 


The fluid-crust limit of this equation does not pose any problem because its derivation does 
not require the lithospheric condition (Eq. ©)• Indeed, the last term of Eq. ( [5^ does 
not diverge in the fluid-crust limit. The fluid-crust limit of Eq. reads 


Ry^iRe) + 6 (n - 1) f e (fe° + 1) = (2n + 1) (l - e (n - 1 - 3f)) . 


(92) 
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5.4 Gravity scaling with crustal parameters 


I recall here the scaling argument explained in Beuthe [2015 , extending it to a crust of finite 
thickness and allowing for a possible density contrast at the crust-ocean boundary. At the 
mantle-ocean boundary {r=Rm), the variables (^ 5 ,^ 7 ) can be related by continuity to the 
six yi solutions within the mantle. In the mantle, the ?/j-vector is a linear combination of 
three independent solutions because there are only three regular solutions at the center of 
the body (if there is a liquid core, the solutions within the mantle can be expressed in terms 


of three unknown constants at the core-mantle boundary [e.g. Sabadini and Vermeersen 
2004| ). 

The three constants of this linear combination reduce to one after applying the free-slip 


condition {y/^{Rm) = 0) and the fluid condition taken in the static limit (Eq. (25)) at the 


mantle-ocean boundary. Both conditions are homogeneous in the sense that they do not 


introduce a constant term that would be independent of the y, (by contrast Eq. (91) is not 
homogeneous). Therefore, the six yi{Rm) at the mantle-ocean boundary and (ys, y^) at any 
radius within the fluid linearly depend on one free constant C, with proportionality factors 
fi{r) depending on the radius and on the structure of the body below the crust (densities, 
radii of interfaces, rheology) but not on the crustal parameters {A,6p°) themselves: 


yi{r) =Cfi{r) 


(93) 


The variables yi{Rm) and (y 5 ,y 7 ) become dependent on {A,5p°) when the constant C is 
determined from the crust-ocean boundary condition on gravity (Eq. (IH))- Now suppose 
that the solutions yi^a = Cafi{r) and yi^ = Cbfi{r) are associated with the crustal parame¬ 
ters {Aa,dp’^) and {Ab,6pl), respectively. Then the ratios yi^a/Uifi are equal (for any i) to 
the ratio ( = Ca ICb- 

Erom the general argument explained above, I can write 


Vb jr) 

2/5(0 


2/7 (?^) 
y?(r) 


for Rm < r < R^ . 


(94) 


The ratio C depends on (A, 5p°) and is related to the yi solutions within the mantle by 


c = 


yi{.Rm) 


(95) 


The ratio C can be related to {kn, k^) by propagating (ys, yg) to the surface. Using Eqs. (89) 

fen + 1 


(90), I write ( at 0(e) as 


C=|l-3e^ 

Pb 


fe ° +1 


(96) 


Eqs. (95) and 


are useful when decomposing tidal heating into crustal and mantle 


contributions (Section |5.6[). 
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5.5 Explicit formulas for kn and hn 

The scaling relation given by Eq. (94) evaluated at the crust-ocean boundary reads 


Vb (Re) ^ Vy (Rs) 
y°,iR,) ~ y°,{Re) 


( 97 ) 


Substituting Eqs. (89)-(92) into Eq. (97), I can solve to 0{e) for kn in terms of /c°. Next, 
I combine this result with the /c„ — hn relation in order to express hn in terms of /i°. The 
resulting formulas read 


fcn + 1 — 


hn — 


1 + 


3$° 

2n+l 


k° + 1 


+ 1) 

hi 


1+1 + 


3g° 


2n+l ri 


/l° ) A + Ay + Hp 


(98) 


where all terms of 0(e) are in the denominator (^ and are interchangeable in these 
terms). The density corrections Kp and Hp are defined below. As the compressibility 
factor does not appear in the formula for kn, crust compressibility has a larger effect 
on hn than on kn (kn weakly depends on crust compressibility through A). 

Density corrections are given by 


Kf, 

Hp 



1 ) {^n + 1 ) “ ( 2 ?^ + 1 ) 



Ap 


2n + l5p 3C 
hi p 2n +1 


(99) 


In the last equation, the first two terms in the right-hand side depend on 5p whereas the 
last term depends on 6p°. 

What is the impact of choosing the ocean density or the physical crust density as the 
fluid-crust density? With the former choice (p° = p), the fluid-crust Love numbers can, for 
example, be computed with the two-layer incompressible body of bulk density made of 
a viscoelastic core and a homogeneous ocean of density p reaching the surface (Eq. (C.3)). 
Beware that it is not correct (if p° = p) to compute the fluid-crust Love numbers with a 
bulk density equal to pi,. Expanding the solution about a zeroth-order configuration of bulk 
density ph affects density corrections of 0(e). This procedure is illustrated in Appendix [Pj 
where I prove that the membrane formulas for Love numbers agree with the analytical 
model of Wahr et al. 2006[ in the rigid mantle limit. 

If the fluid-crust density is equal to the physical crust density (p° = p), then 6p° = 0 
and Kp = 0. The resulting formula for kn is formally identical to the one derived in 
the massless approach (Eq. (57) of Beuthe 2015| ). Density corrections, however, are now 
hidden in A:°: when computing A;°, one must take into account an ocean of density p and 
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a fluid-like crust of density p. For example, Eq. (C.9) gives the fluid-crust Love numbers 


for the three-layer incompressible body made of an infinitely rigid mantle, a homogeneous 
ocean, and a homogeneous fluid crust differing in density. 


5.6 Micro-macro equivalence in tidal dissipation 

The imaginary part of Love numbers is much smaller than their real part and is thus 
difficult to measure with geodetic methods (altimetry or gravity). The imaginary part 
of /c 2 , however, manifests itself in a very different way: the tidal energy dissipated by 
viscoelastic friction in the whole body is proportional to Im{k 2 )- Consider a synchronously 
rotating body with spin rate cu, orbital eccentricity e and obliquity I. The global heat flow 
due to tidal dissipation is given by 


■ 5 [ojRY’ ^ /.NT 


( 100 ) 


where Tg = (21/5)e"^-|-(3/5) sim / (Eqs. (41)-(42) of Beuthe 2013 ). This macro approach 


to tidal dissipation is equivalent to the micro approach, in which the dissipated energy is 
computed by integrating over the whole body the product of the microscopic stress and 
strain rate. In particular, the global heat flow due to dissipation in a thin floating crust is 
given by 

3 (ioR)^ 

2 G 


E. 


crust — 


C|h2p/m(A) Tg. 


( 101 ) 


This formula was derived in the massless membrane approach (Eq. (98) of Beuthe\ [2015| ) 
but one can prove that it still holds in the massive membrane approach. 

If the core and mantle are purely elastic, the global heat flow shoul d be equal to the 


2015 


to 


heat flow coming from the crust: E = Ecrust- Applying the method of \Beuthe 
Eq. (98), I decompose Im(k 2 ) into contributions from the crust and core-mantle system 
(denoted c-m)\ 

Im{k2) = [Im{k2)]^^^,^ + [Im{k2\_^ , ( 102 ) 


where 


VMk2)]crust = \h2? Im{K) ^ 

[Im{k 2 )]c-m = , 


(103) 


in which is defined by Eq. ( |95| ). That [Im{k 2 )]crust is indeed the crustal contribution 
can be checked by substitution into Eq. ( |100 ): the result is Eq. (101) as it should be. 
With respect to the fluid-crust model, the core-mantle contribution of the physical model 
is reduced by the factor where C, is the reduction in radial displacement of the mantle- 
ocean boundary due to the viscoelasticity and to the density contrast 5p° of the crust. 
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6 Numerical benchmarking 


In membrane formulas, there is a clear distinction between the contributions of the shallow 
interior (crust and crust-ocean boundary) and of the deep interior (the rest). The latter 
influences Love numbers through the fluid-crust Love numbers /i° and The accuracy 
of the membrane formulas thus depends on two different things: 

1. the crust thickness. This type of error is intrinsic (membrane approach = perturbative 
expansion) and unavoidable (higher-order corrections due to finite crust thickness). 

2. the choice of a particular fluid-crust model. This type of error can be reduced to 
zero; it depends on the degree of complication one is willing to accept in order to 
compute fluid-crust Love numbers. 


In this section, I quantify the intrinsic error of the membrane formulas in specific models 
of the shallow interior of Europa and Titan, postponing the error analysis of fluid-crust 
models to Sectionj^ Benchmarking is done with the program love.f included in the software 
SatStress (available at http://code.google.eom/p/satstress) [ Wahr et al. I also use 

the propagator matrix method in the Fourier domain. 

As a first step, I benchmark the membrane formulas for tidal Love numbers (radial and 
gravitational). Titan is a good choice for this purpose because its Love number k 2 was re¬ 
cently estimated with Cassini data \Iess et ^ 2012 ; this is actually the only measurement 
of a tidal Love number for an icy satellite. As a second step, I will benchmark the mem¬ 
brane formula for the tilt factor. Europa is chosen here for two reasons: (1) the tilt factor 
can be used to estimate the crust thickness [ Wahr et al. 2006 which is a key parameter 
regarding Europa’s habitability, and (2) convection may occur in Europa’s crust, making 
it a good laboratory for the analysis of the influence of crustal rheology on Love numbers. 


6.1 Models of Europa and Titan: shallow interior 

6.1.1 Density of crust and ocean 

Testing the intrinsic error of membrane formulas does not call for detailed models of the 
deep interior: core and mantle are treated together as an incompressible homogeneous 
layer with fixed radius {core-mantle system) while the ocean is approximated as a layer of 
uniform density. The parameters describing the shallow interior are the crust-ocean density 
contrast and the thickness, density, and rheology of the crust. For each set of shallow 
interior parameters, I adjust the density of the core-mantle system so that the bulk density 
is equal to the physical value (Table [^ . I do not try to fit the moment of inertia because 
it does not lead to more realistic models as long as the core is not distinguished from the 
mantle. 

What do we know about the density of the ocean? An ocean made of pure water has 
a density of 1000 kg/m^ at 273 K and atmospheric pressure. Pressure effects increase the 
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density by about 0.45 kg/km for Titan (e.g. Fig. 1 of Mitri et al. [2014] ), with a comparable 
effect on Europa due to the similar surface gravity. For Titan, this means that the mean 
density of a 400km-thick pure water ocean below a 50km-thick crust is close to 1100 kg/m^. 
For Europa, the ocean is not as deep and the mean density of a pure water ocean is less than 
1050 kg/m^. Besides the pressure effect, the ocean density is affected by the presence of 

; ammonia can lower the uncompressed density down to 950 kg/m^ 


solutes Sohl et al. 


2010 


whereas dissolved salts increase the uncompressed density up to 1200 kg/m'^. Fortes 2012 


uses these bounds to construct two models of Titan in which the compressed mean density 
of the ocean is either 1020 kg/m^ (‘light-ocean’) or 1280 kg/m^ (‘dense-ocean’), assuming 
a lOOkm-thick crust and a 250km-thick ocean. Mitri et al. [2014| obtain a similar upper 
bound by imposing that the density at the bottom of Titan’s ocean is less than the density 
of the high-pressure ice layer below. While the compression effect is smaller for Europa, 
Kargel et al. [2000 consider even higher density solutions for Europa’s ocean. It is thus a 
reasonable choice to adopt 1020 and 1280 kg/m^ as lower and upper bounds for the ocean 
density in Europa and Titan. 

What do we know about the density of the crust? Pure water ice at atmospheric pres¬ 


sure has a density varying between 917kg/m^ at 273 K and 933kg/m^ at 100 K Feistel 
and Wagner, 2006 , providing a lower bound on the crust density (set here at 930 kg/m^). 


Porosity is not taken into account because it probably only affects a thin layer near the 
surface. The crust, however, is probably highly impure and chemically layered: ice in¬ 


cluding silicate dust or salt hydrates could have a density of 1050 kg/m^ Schubert et al. 
2009|, while Kargel et al.\ [2000 and Spaun and Head 


2001 consider eutectic mixtures 


with magnesium sulfates having a density of 1144kg/mT Since the membrane formulas 
depend on the crust density through the density contrast 5p/p., I adopt a slightly higher 
upper bound of 1167 kg/m^ for the crust density. With this choice, the density contrast 
5p/p is the same whether the crust and ocean densities take the lower values (930 and 
1020kg/m^) or the higher values (1167 and 1280kg/m^). Table summarizes the three 
test cases (‘Light’, ‘Mixed’, and ‘Dense’). 

Table 6: Density models for the crust and ocean of Europa and Titan. The density of the core¬ 
mantle system is adjusted so that the bulk density remains the same. The top of the mantle is at 
a depth of 170 km for Europa and 350 km for Titan. 


Density model 

Crust 
p (kg/m3) 

Ocean 
p (kg/m^) 

Contrast 

dp/p 

Light 

930 

1020 

-0.09 

Mixed 

930 

1280 

-0.27 

Dense 

1167 

1280 

-0.09 
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6.1.2 Rheology 

Table gives the values adopted for the viscoelastic parameters in the core-mantle system 
(purely elastic and incompressible) and in the crust. The value chosen for the shear modulus 
of the core-mantle system is justified in Section The critical viscosity r]crit is defined by 
Eq. (A.4). The rheology of the crust is discussed in more detail below. Although ocean 


compressibility is listed as an input of SatStress, Saito 1974 showed that the stress-strain 


relation of a liquid layer is irrelevant when computing static deformations. Thus SatStress 
output is in principle independent of ocean compressibility which can be set either to its 
physical value (about 2 GPa) or to a much larger value simulating incompressibility. 

Table 7: Viscoelastic parameters (viscosity is only relevant to Europa; Titan is treated as an 
elastic body). 


Parameter 

Symbol 

Value 

Unit 

Shear modulus of core-mantle system 

l^m 

40 

GPa 

Bulk modulus of core-mantle system (for SatStress) 

Km 

102O 

Pa 

Shear modulus of elastic ice 

g-E 

3.5 

GPa 

Poisson’s ratio of elastic ice 

ee 

0.33 

- 

Viscosity of top ice layer“ 

‘Htop 

10" 

^crit 

Viscosity of bottom ice layer® 

‘Hbot 

10/1/0.1 

^crit 


Europa: r/crit = 1-71 X 10^'* Pa.s; Titan: rjcrit = 7.67 X 10^^ Pa.s (Eq. i A.4 i). 


Viscous effects are probably small within Titan’s crust. Uncompensated gravity anoma¬ 


lies indeed suggest that the crust is not convecting Hemingway et al ., 2013; Mitri et al. 


2014 Lefevre et al. 2014| . I will thus assume that Titan’s crust is purely elastic. There 
are no comparable gravity data constraining the state of Europa’s crust which could be 
either in a conductive regime or in a stagnant lid regime (Section 3.4). In Moore 2006[ ’s 
study, Europa’s crust always convects and the total crust thickness d varies between 20 km 
and 120 km, depending on the viscosity of the convecting layer. The three following 
cases can be distinguished (recall that the critical state is defined by the critical viscosity 
Vcrit = see Table 0 and Appendix [A]) . If the convecting layer is close to its critical 

state, the crust is thin [d ~ 20 km) and the top ice layer is relatively thick (say 2/5 of d). 
If the convecting layer is elastic-like {rjbot dcrit), the crust can be thick {d up to 100 km) 
and the thickness of the top ice layer scales with the total thickness (it is thus relatively 
thick as in the critical case). If the convecting layer is fluid-like {r]bot 'C rjcrit)-, the crust 
can be thick {d up to 120 km) but the top ice layer does not scale with d: its thickness 
remains close to its value at critical viscosity. 

Table summarizes the three test cases for crustal rheology. The values of the effective 
viscoelastic parameters {fL,u,x) are computed with Eqs. (42), (43) and (45) using Maxwell 


rheology (Eq. (A.2)). In the elastic-like case, the effective parameters are nearly equal to 
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their elastic value. In the critical model, the shear moduli of the conductive and convective 
layers are equal to /ig and + ^)/2, respectively, so that their weighted average is 

0.7+0.3h In the fluid-like case, the real part of the effective shear modulus is approximately 
equal to the elastic shear modulus multiplied by the relative thickness of the lithosphere 
(Eq. (49)), while the effective Poisson’s ratio is getting closer to its value in the fluid limit 
{v = l/2). 


Table 8: Rheology models for Europa’s crust. The crust (of thickness d) is made of two homoge¬ 
neous layers: the top layer (of thickness dtop) is purely elastic while the bottom layer is viscoelastic. 


Rheology model 

'Hbot / ^crit 

dtop / d 

fi/dE 

u 

X 

Elastic-like 

10 

0.4 

0.994-t 0.059 i 

0.331 -0.009 i 

0.506-t 0.020* 

Critical 

1 

0.4 

0.700-t 0.300 i 

0.385-0.034 i 

0.415-t 0.140* 

Eluid-like 

0.1 

0.1 

0.109-t 0.089* 

0.435-t 0.061* 

0.061-t 0.067* 


6.2 Love numbers of Titan 


In this section, I compute the gravitational and radial Love numbers for models of Titan 
differing in ocean and crust densities (Table . The ocean density is a crucial parameter 
which could account for the large measured value of k 2 \Iess et ^ [2012] . I assume that 
the crust is purely elastic (Section |6.1.2D. 

Fig.0 shows the Love numbers /c 2 and /12 of Titan in terms of the relative crust thick¬ 
ness. Solid curves show the membrane predictions (Eq. (98) with an incompressible purely 


elastic core; the fluid-crust model is given by Eq. ( |C.3 )). Dashed curves show SatStress 
predictions if the core-mantle system is incompressible. Dotted curves show the results 
of the propagator matrix method (equivalent to SatStress predictions for an incompress¬ 
ible body). Membrane and SatStress predictions agree perfectly in the limit of zero crust 
thickness. For non-zero crust thickness, the agreement is excellent up to a relative crust 
thickness of 5% and reasonable up to 10%. Compared to the membrane approach, the 
accuracy of the propagator matrix method is similar for k 2 but much worse for /i 2 , because 
crust compressibility has a much larger on /12 (up to several percents) than on k 2 - 

Fig. [6] shows the relative error between the membrane and SatStress predictions. Solid 
curves show the error for the Love numbers of Fig. (for which the core-mantle system is 
incompressible). The error is less than 1% (for ^ 2 ) and 0.5% (for /i 2 ) if the relative crust 
thickness is less than 5%. 


6.3 Tilt factor of Europa 


In Section 6.2,1 showed that membrane formulas approximate very well k 2 and / 12 , and that 
the accuracy improves as the crust thickness d decreases. These observations, however, do 
not prove that the membrane formulas are the correct perturbative expansions of thick shell 
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Figure 5: Tidal Love numbers of Titan as functions of relative crust thickness: (A) k 2 , and (B) 
/i 2 . Solid curves are membrane estimates. Dashed curves are SatStress predictions. Dotted curves 
are predictions made with the propagator matrix method. Light/Mixed/Dense refer to the density 
models of Table See Section 


6.2 


for details. 


A B 




Figure 6: Tidal Love numbers of Titan: relative error of the membrane prediction for (A) k 2 , and 
(B) /i 2 . The error corresponds to the difference between the solid and dashed curves of Fig. 
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theory: it only shows that all crustal contributions to Love numbers are proportional to d. 
The best way to test the dependence of Love numbers on crustal parameters is to study the 
tilt factor 72, because this quantity vanishes in the limit d —)• 0. In particular, 72(0) (the 
slope of the tilt factor at d = 0) is a good indicator of the correctness and completeness of 
membrane formulas. For example, Beuthe 2015 concluded from benchmarking 72(0) that 


compressibility terms are missing in the k 2 — ^2 relation derived from classical thin shell 
theory. 

In the membrane approach, 72 (0) is given by the right-hand side of Eq. (69) divided by 
the relative crust thickness e. Fig. shows the real and imaginary parts of 72(0) for the 
critical rheology model of Europa (Table [^. The ticks on the x-axis correspond to four 
different approximations of the tilt factor (denoted by circles, squares and triangles): 


1 . 

2 . 

3. 

4. 


A/i2: dominant elastic term (thin shell theory prediction of Beuthe\ [2015]), 


(A -|- A^)h 2 - full elastic term, including the compressibility correction. 


A7’/i2: terms proportional to /12, including the minor density correction Ap, 

72: full membrane prediction, including the major density correction —5{6p/p). 


Fig.0 demonstrates that all contributions are required in order to reach a near perfect 
agreement with the predictions of the numerical benchmark (horizontal dashed lines). Re¬ 
garding the real part of 72(0), one can say that 


• the slope of the tilt factor decreases as the crust becomes softer (elastic-like —)• critical 
—)■ fluid-like); 

• the elastic term is generally dominant though the density term can become large if 
the lower crust is fluid-like, at least for large satellites (pE/ipgR) ~ 1); 

• compressibility and minor density corrections are negative while the major density 
correction is positive: the various corrections thus partially cancel each other; 

• if the rheology is elastic-like, the tilt factor is not very sensitive to the ocean density 
(Light and Mixed models yield similar results); this is also true of Titan, but not of 
all icy satellites because it depends on the values of the crust parameters and of the 
surface gravity; 

• if the rheology is fluid-like, the tilt factor is sensitive to the crust-ocean density 
contrast, but does not depend much on the crust and ocean densities taken separately. 


Regarding the imaginary part of 72(0), it is significant (with respect to the real part) if 
the rheology is critical or fluid-like. The compressibility correction has a 10% effect while 
density corrections contribute little to /m(72(0)). 
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r2'(0) Re 72'(0) 


ELASTIC-LIKE 


CRITICAL 


FLUID-LIKE 



A /?2 (A+A^)/72 /^jh2 y2 ^ ^2 (A+A^)/72 A7'/72 y2 ^ ^2 (A+A^)/72 l\.jh2 

Various approximations of tilt factor Various approximations of tilt factor Various approximations of tilt factor 


Figure 7: Slope (at zero crust thickness) of Europa’s degree-two tilt factor: Re(72(0)) (upper pan¬ 
els), and Im(72(0)) (lower panels). The different columns correspond to the crustal rheology models 
of TableCircles/squares/triangles are membrane estimates for the Light/Mixed/Dense models 
of Table ^ Empty symbols show partial membrane estimates (with or without compressibility and 
density corrections). Filled symbols show the full membrane estimates (Eq. (69)). Dashed lines 
show SatStress predictions (d = 1km), with L/M/D standing for Light/Mixed/Dense models. Eor 
membrane estimates, h 2 is evaluated with Eq. fcTt : h 2 = 1.27 (Model L) or h 2 = 1.35 (Models M 


and D). See Section 6.3 for details. 
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Consider now the tilt factor at non-zero values of the crust thickness. Theoretically, 


the right-hand side of Eq. (69) could be evaluated with /12 equal to its fluid-crust value (of 
0(1)) because it is already multiplied by a term of 0(e). However, the 0(e) terms in /i2 
have a large multiplying factor so that /12 is rather sensitive to crustal thickness (see Fig. 
for Titan). By contrast, next-to-leading corrections to A are relatively much smaller. I 
can quantify this assertion with the homogeneous crust model of Appendix For static 
tides, 72 = A/i2 and /12 = + ^2r^) where A = Xafi and Xa ~ (24e/ll)(l -|- 4e/ll) 


(see Table 11). At next-to-leading order, A is modified by the factor 1 -|- 4/lle ~ 1 -|- 0.4e 
whereas /12 is modified by the factor l + N^^k ~ l-|-4.7e (assuming /i2^ ~ 1.25 and jl ~ 1.71 
for Europa). If Europa’s crust is 40 km thick, next-to-leading order corrections modify A 
and /i2 by 1% and 12%, respectively. Therefore the membrane estimate of the tilt factor 


is much more accurate if /12 is evaluated at 0(e) in the right-hand side of Eq. (69). 

Fig. shows the tilt factor of Europa as a function of the relative crust thickness for 
the density and rheology models of Tables and According to the discussion above, the 
right-hand side of Eq. (69) is evaluated with /i2 depending on the crust thickness (Eq. (98)). 
This effect accounts for the concavity of the tilt factor curves in Fig.[^ Membrane estimates 
(solid curves) are in good agreement with SatStress predictions (dashed curves). As noted 
for Fig. the tilt factor becomes smaller if the lower part of the crust becomes fluid-like. 

Fig.j^shows the error for the various density and rheology models. The error is smallest 
(resp. largest) if the rheology is elastic-like (resp. fluid-like) and the crust-ocean density 
contrast is small (resp. large), because next-to-leading corrections are larger for density 
terms than for elastic terms. The error on the tilt factor is one order of magnitude larger 
than the error on Love numbers because the tilt factor is a small quantity (of 0{e)) resulting 
from the cancellation of the dominant terms (of 0{1)) in the Love numbers. 


7 Influence of deep interior 

7.1 Fluid-crust models of Europa and Titan 

In the membrane approach, the deep interior affects the Love numbers through the fluid- 
crust Love numbers /i° and fe°. Fluid-crust models are thus sufficient when studying the 
influence of the deep interior. It is easy to add the effect of the crust by using the membrane 
formulas (Eq. ([9^). In the fluid-crust limit, large icy satellites can be roughly described 
with four layers (not necessarily all present): high-density core, silicate mantle, layer of 
high-pressure ice, and surface ocean. There are however only two observational constraints 
on the interior structure: the total mass and the moment of inertia (Table [^. 

For Europa, I consider three-layer models made of an iron core (solid or liquid, light 
or dense), a silicate mantle and a pure water ocean. Similar models have been used by 
Anderson et al. [1998 , Moore and Schubert 2000] , Tobie et al. [2005 , and Schubert et al. 

2009 . The density of the three layers are fixed as in Table and the radii of the interfaces 


are computed so that the constraints of total mass and moment of inertia are satisfied (this 
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Figure 8: Europa’s degree-two tilt factor (absolute value) as a function of crust thickness. Solid 
curves are the membrane estimates (Eq. ([6^ in which /12 is evaluated with Eq. (98)). Dashed 


curves are SatStress predictions. The different panels correspond to the crustal rheology models 
of Table L/M/D stand for Light/Mixed/Dense denoting the density models of Table See 


Section 6.3 for details. 



Figure 9: Relative error (absolute value) on Europa’s degree-two tilt factor in the membrane ap¬ 
proach. The error is defined as the relative difference between membrane and SatStress predictions. 
These errors are a bit larger than those deduced from Fig. because they include the imaginary 
part. Solid, dotted, and dash-dotted curves correspond to the elastic-like (E), critical (C), and 
fluid-like (F) cases, respectively. L/M/D identify the density models as in Fig. The error curves 
EL and ED nearly coincide. 
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makes sense if the solid crust to be added to the fluid-crust model has the same density as 
the ocean). When studying the influence of the liquid core, I will also consider a range of 
densities for the silicate mantle. 

For Titan, I consider three-layer models made of a silicate core, a mantle of high- 
pressure ice, and an ocean. As the total mass and moment of inertia are not very sensitive 
to the thickness of the high-pressure ice layer, I construct fluid-crust models on the basis 
of the interior models of Mitri et al. [2014] (pure water ocean and high-density ocean, both 
with a 50 km-thick crust, see their Fig. 1). I approximate the core, mantle and ocean as 
homogeneous layers and I replace the crust by an ocean layer (the total mass and moment 
of inertia thus slightly differ from the observed values). The density structure of the two 


models is given in Table 10 


Regarding elastic parameters, all layers are incompressible so that only the shear modu¬ 
lus must be specified. The shear modulus of Europa’s iron core is zero if it is liquid, 75 GPa 
if it is made of solid FeS (low density), and 100 GPa if it is made of solid Fe, respectively 
Sohl et al., 2002] . The silicate mantle has a shear modulus between 40 GPa (hydrated sili¬ 


cates of low density) and 80 GPa (olivine). High-pressure ice has a shear modulus between 
6 GPa (Ice V) and 8 GPa (ice VI) \Sotin 1998 


Tidal deformations are computed with the propagator matrix method. The resulting 
formulas for the three-layer models are complicated and will not be given here. I will show 
that it is a good approximation to consider the core and mantle together as a layer of 


uniform density, so that one can resort to much simpler two-layer formulas (Eq. (G.3)) 
For this purpose, the mean density of the core-mantle system is given in Tables and 10 


The viscoelastic effect of the deep interior being rather small, it is interesting to define 
the change of the Love numbers with respect to the rigid mantle model, i.e. the quantities 


' k° — h° h° — h' 

[ek eh) = ' 


k: 


2r 


k 


2r 


(104) 


where (fe2r)^2r) given by Eq. (C.l). Since /12 = -b 1 and = (3^°/5)/i2r, and eh 
are related by 

Cfc = ^ Eh • (105) 

It is thus sufficient to compute either or eh- For Europa, e^, ~ 5e/i whereas ~ (5/2)e/i 
for Titan. 


7.2 How good is the rigid mantle approximation? 

At several occasions in this paper, I use the rigid mantle approximation, i.e. the limit of 
infinite mantle rigidity. Though the mantle in icy satellites is not more rigid than Earth’s 
mantle in an absolute sense, it is very rigid with respect to the ocean layer that surrounds 
it. The largest part of the tidal response occurs in the ocean layer and in the crust. I 
will quantify this statement in two ways, by computing the effect of an elastic mantle on 
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Table 9: Fluid-crust models of Europa: density structure. Crust and ocean are assumed to have 
the same density so that the bulk density and the moment of inertia of the fluid-crust model are 
equal to the observed values. 



Light 

core 

Dense core 


radius 

density 

radius 

density 


(km) 

(kg/m3) 

(km) 

(kg/m^) 

Core 

727 

5150 

518 

8000 

Mantle 

1445 

3300 

1449 

3300 

Core-mantle system 

1445 

3535 

1449 

3515 

Ocean 

1560.8 

1000 

1560.8 

1000 

Table 10: Fluid-crust models of Titan: 

density structure. 

The two models result from the homo- 

geneous layer approximation of the end member models of Mitri et al. 

2014 . Since the ocean layer 

that replaces the crust {d = 50 km, pc = 

935kg/m^) does not have the same density as the crust. 

the bulk density and the moment of inertia (not given here) slightly differ from the observed values. 


Light 

ocean 

Dense ocean 


radius 

density 

radius 

density 


(km) 

(kg/m3) 

(km) 

(kg/m^) 

Core 

2094 

2542 

1968 

2645 

Mantle 

2134 

1367 

2183 

1373 

Core-mantle system 

2134 

2477 

2183 

2305 

Ocean 

2574.8 

1118 

2574.8 

1270 

Whole body 

2574.8 

1892 

2574.8 

1901 

the deformations of the mantle and crust. As 

an illustration, I choose the solid light-core 


the shear modulus is the same in the core and mantle. For Titan, the shear modulus in 
the rocky core and in the high-pressure ice layer differ by a factor of ten. 

First, I compute the relative deformation of the top of the mantle with respect to the 


surface, i.e. the ratio /h^ (given for the two-layer model by Eq. (C.8)). Fig. 104 shows 
how varies with the core rigidity. Solid (resp. dashed) curves show results for the 

two-layer (resp. three-layer) model. Plausible values for the shear modulus of the core 
range from 40 to 70GPa (shaded band) so that h^/h^ is about 1% for Titan and 2% for 
Europa. The deformation of the mantle becomes large (say h^/h^ > 20%) if the shear 
modulus of the core decreases below the value for ice I (between 2 and 4GPa). Below a 
shear modulus of about 0.2 GPa, h^/h^ tends to a constant value corresponding to the 
deformation of a fluid body. 
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Second, I compute the change of the surface deformation with respect to the rigid 
mantle model, i.e. the quantity eh given by Eq. (104). Fig. [lo^ shows how this quantity 
varies with the shear modulus of the core. Conventions are the same as in Fig. [TgA: 
solid (resp. dashed) curves show results for the two-layer (resp. three-layer) model. For 
plausible values of the core rigidity (shaded band), the elasticity of the mantle increases 
the surface deformation by 0.2-0.4% (for Titan) and by 0.7-1.2% (for Europa). Similarly 
to what was observed for mantle deformation, the elasticity of the mantle starts to have a 
large effect when the shear modulus of the core decreases below the value for ice I. Below 
a shear modulus of about 0.2 GPa, eh tends to a constant because tends to its fluid 
Love number value. Fig. [T0) G shows the error made by using the two-layer model instead 
of the three-layer model. For plausible values of the core rigidity (shaded band), the error 
made by treating the core and mantle as one homogeneous layer is one order of magnitude 
smaller than the correction due to the elasticity of the core and mantle. 

It thus makes sense to consider the core of a large icy satellite as being infinitely rigid 
if its rigidity is larger than ~ 4 GPa, as fluid if its rigidity is less than ~ 0.2 GPa, and as 
elastic in-between. For smaller icy satellites, the approximation of an infinitely rigid core 
is even better because Love numbers depend on the nondimensional shear modulus of the 
mantle which is inversely proportional to the product of the surface radius and surface 


gravity (Eqs. (G.3)-(G.4)). The other lesson of Fig. is that the two- and three-layer 


models give very similar results, especially for large (and thus plausible) values of the core 
rigidity. The density contrast between core and mantle is thus a secondary factor, as is the 
presence within Titan of a high-pressure ice layer with a low shear modulus (the layer is 
too thin to affect much the response of the core). 


7.3 Influence of the liquid core 

If Europa’s iron core is liquid, how much will it increase the surface deformation? For the 
model with a solid core, I concluded that the density contrast between core and mantle is 


a secondary effect (compare solid and dashed curves in Fig. 10). This observation suggests 
the following approximation: neglect density effects by assuming that core and mantle 
have the same density pcm- In that case, the core-mantle boundary does not interact 
gravitationally with the ocean, so that one can treat the core and mantle as one layer 
(core-mantle system) with an equivalent shear modulus denoted pcm- 

How is Pcm related to pm, the shear modulus of the mantle? Let us first ignore the 
ocean. The liquid core-solid mantle system is a two-layer body of density pcm, surface 
gravity pm and radius Rm- The radial Love number of this model is given by Eq. (E.8): 




1 


2 1 -|- 2 pm 


(106) 


where pm = Pm/{PcmdmRm) and Xa is a function of x = Rc/Rm defined by Table 11 
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Figure 10: Influence of core-mantle elasticity on tidal deformations of Europa and Titan: (A) 
relative deformation of the top of the mantle with respect to the surface, (B) change of the surface 
deformation with respect to the rigid mantle model (Eq. (104)), and (C) error on /12 between two- 
and three-layer models. Thick curves refer to Europa (light-core model) while thin curves refer to 
Titan (dense-ocean model). Solid (resp. dashed) curves show results for the two-layer (resp. three- 
layer) model. All models have a fluid crust. Shaded bands indicate plausible ranges for the core 
rigidity. Dotted vertical lines separate fluid, elastic, and rigid regimes for the core. In panel (C), 
the convergence of the curves in the fluid limit is a coincidence. See Section 7.2 for details. 


Equating Eq. (106) with the Kelvin-Love formula, /12 = (5/2)/(l -|- (19/2)/ic 


see 


Eq. (E.9)), I define the equivalent shear modulus of the core-mantle system by 


^^cm — l^'m ■ 


(107) 


Fig. [TI] A shows the ratio ^cm/lJ-m in terms of the relative size of the core with respect 
to the mantle (Rc/Rm)- It varies between one (no liquid core) and zero (no mantle). If 
the core extends to about one-half the radius of the mantle, ficm ~ Q-hhUrn so that the 
shear modulus of the mantle (70GPa) is reduced to an equivalent shear modulus of about 
40 GPa. 

The concept of equivalent shear modulus remains valid if there is an ocean above the 
mantle, as long as the core and the mantle have the same density (this can be proven with 
the propagator matrix method). Fig. shows the change of the surface deformation 
with respect to the rigid mantle model (see Eq. (104)) as a function of the relative size of 
the liquid core. The solid curve corresponds to the results of the two-layer model made 
of a viscoelastic mantle below an ocean (Eq. ( C.3[ )) in which the shear modulus of the 
mantle (70 GPa) is replaced by the equivalent shear modulus of the core-mantle system. 
Big dots show the results for the three-layer models of Table with the shear modulus of the 
mantle equal to 70 GPa. Dashed curves show the results for three-layer models satisfying 
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the constraints of total mass and moment of inertia, in which the mantle density varies 
between 2500kg/m^ (large core) and 3800kg/m^ (small core); the core and ocean densities 
are the same as in Table (these models are inspired by the models of Table 3 in \SchubeH 
2009] ). If there is no liquid core, eh ~ 0.7% (as in Fig. 


et al. 2009] ). If there is no liquid core, eh ~ 0.7% (as in Fig. jlOp with //core = 70GPa). 
The presence of a liquid core increases the surface deformation up to eh ~ 1.7% if the core 
and the mantle are both light. This figure also shows that the results of three-layer models 
do not differ much from those of the two-layer model using the equivalent shear modulus. 
As in the solid core case, the density stratification of the core-mantle system is a secondary 
factor. 


A B 




^c/ ^c/ 


Figure 11: Influence of a liquid core on the surface deformation of Europa: (A) equivalent shear 
modulus of the core-mantle system, and (B) change in surface deformation. The a;-axis variable is 
the relative size of the liquid core with respect to the mantle. In panel (A), the y-axis variable is 
the equivalent shear modulus of the core-mantle system (Eq. (107)) divided by the mantle shear 
modulus. In panel (B), the y-axis variable is the change in surface deformation with respect to the 
rigid mantle model (Eq. (104)). The solid curve shows the two-layer model made of a viscoelastic 
mantle below an ocean (Eq. (C.3)) in which the shear modulus of the mantle is replaced by the 
equivalent shear modulus of the core-mantle system. Big dots show the results for the three-layer 
models of Table ]^ while dashed curves show the results of three-layer models for a wide range of 
mantle densities. All models have a fluid crust. See Section 


7.3 


for details. 


7.4 Ocean stratification and screening effect 


For a satellite with a subsurface ocean, the major parameter determining the magnitude 
of Love numbers is the ocean-to-bulk density ratio ^ (or in the corresponding fluid- 
crust model). This is easily seen with the two-layer body made of a rigid mantle and an 


incompressible homogeneous ocean, for which Eq. (C.l) gives 


k° — 

'^2r — 


3^ 

5 - 3 ^° 


(108) 
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Figure 12: Screening effect due to ocean density stratification. The model is a three-layer body 
made of an infinitely rigid mantle M below two ocean layers L and U. In panel (A), the ocean layers 
have the same density and the rigid mantle acts as a gravitational brake on the surface deformation. 
In panel (B), the lower ocean layer has the same density as the mantle and screens its gravitational 
braking effect. 


This formula yields values in the range [0.25, 0.35] for Europa and in the range [0.47, 0.71] 
for Titan (assuming an ocean density between 1000 and 1300kg/m^). 

What happens if the ocean is stratified in density? For a hydrostatic body, increasing 
density concentration always decreases (fluid) Love numbers; this mass concentration effect 
can be proven with the Radau relationship which relates the gravitational fluid Love number 
to the moment of inertia [e.g. Schubert et al., 2004| . By contrast, Mitri et al. [2014 observe 
that density stratification (mainly due to pressure) within Titan’s ocean leads to a 3 to 4% 
increase in k 2 - Therefore density stratification has different effects in hydrostatic bodies 
and in bodies with elastic layers. The difference can be attributed to a screening effect, 
which I will first explain qualitatively before proceeding to numerical estimates. 

Consider the two-layer body made of a viscoelastic mantle and a surface ocean whose 
Love numbers are given in Appendix 0 One expects that the surface deformation gets 
smaller if the mantle rigidity gets larger, and this can be indeed proven with Eq. (C.3). 
The reason is that a more rigid mantle deforms less and thus contributes less to the induced 
geoid perturbation: it acts as a ‘gravitational brake’ on the deformation of the ocean layer. 
The effect vanishes if the mantle and the ocean have the same density, in which case the 
mantle has no effect whatsoever on the surface deformation. This observation suggests a 
way to cancel the gravitational braking effect by redistributing mass from the top to the 
bottom of the ocean. For example, insert at the mantle-ocean boundary a thin liquid layer 
having the same density as the mantle; decrease also slightly the ocean density so that 
the total mass remains constant. The density distribution is nearly unchanged so that 
the mass concentration effect is not significant. The deformation of the thin liquid layer 
completely screens the gravitational braking of the mantle, effectively transforming it into 
a fluid layer (Fig. 12). As a result, the surface deformation increases and so do the Love 
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numbers. In more realistic models (Section |7.1| ), the core-mantle system is not of uniform 
density and the density at the bottom of the ocean cannot be larger than the mantle 
density. Therefore the screening effect is partial: what matters is the density contrast 
between the mean density of the core-mantle system and the bottom layer of the ocean. 

The screening effect can be quantihed with the three-layer body made of an infinitely 
rigid layer M (for mantle) and two ocean layers L (for lower) and U (for upper). Layer 
densities are denoted {pm-, Pl, Pu) while layer radii are denoted {Rm, Rl,Ru)- This model 
is equivalent to the three-layer fluid-crust model of Appendix [C| in which the upper ocean 
layer takes the place of the fluid crust (Eq. (C.9)). If the densities of the two ocean layers 


are equal, the Love numbers reduce to those of a two-layer body with a rigid mantle 


(Eq. (108)). I parameterize the density stratification as follows: 


PL 

PU 


= p + 6p, 

= p — aSp . 


(109) 


where <5/? > 0 and a > 0 are free parameters. Imposing that p is the average density of the 
whole ocean, I obtain a constraint on the boundary between the upper and lower ocean 
layers: 

1/3 

, ( 110 ) 


(xM + a\ 


where xl = Rl/Ru and xm = Rm/Ru- If a 1, the boundary is located close to the 
mantle while a = 1 corresponds to ocean layers of equal volume. 

I consider two stratification models: the first one is an nearly linear increase in density 


with depth through the whole ocean, due to a pressure effect, as in Mitri et al. 2014 


The value a = 1 is a staircase approximation of this variation. Another kind of ocean 
stratihcation arises from the presence of a thin layer of higher density at the bottom of the 
ocean. This case is modeled with the value a = 0.01. 

The Love number can be computed with Eq. (C.9) in which {x,^°,^°) are replaced 
by {xL,PL/Pb,Pu/Pb) with {xl,Pl,Pu) given by Eqs. (|109|)-([ITol). Eig. shows the 


dependence of on the density deviation 6p for the two Titan models of Table [lb| (core 
and mantle are treated as one rigid layer). Solid (resp. dashed) curves correspond to the 
case a = 1 (resp. a = 0.01). The starting values of k^ (at 6p = 0) differ greatly between the 
two models because the mean density of the ocean is different. As the ocean becomes more 
stratihed, k^ increases with an nearly linear dependence on Sp. Eor these models, the slope 
of the curve does not depend much on the thickness of the lower layer (i.e. the parameter 
a). There is however an important difference between the cases a = 1 and a = 0.01. In 
the former case, 5p is about 80 — 85kg/m^ (see Eig. 1 of Mitri et al. 2014 ) and cannot 


be much larger otherwise either the density at the top of the ocean is lower than the crust 
density, or the density of the ocean bottom is higher than the mantle density. This kind of 
density stratification thus increases k"^ by no more than 3%. In the latter case (a = 0.01, 
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Figure 13: Influence of ocean density stratification on the gravitational Love number: (A) as a 
function of the density deviation of the lower ocean layer, and (B) change of as a function of the 
mantle radius if 6p = lOOkg/m^. The Love number k^ is computed with a three-layer model made 
of a rigid mantle and two ocean layers (Eq. (C.9)). The lower ocean layer is either thick {a = 1, 


solid curves) or thin (a = 0.01, dashed curves). Panel (A) shows k^ for the two models of Table 10 


Panel (B) shows the change in k^ for a suite of models parameterized by their mantle radius (the 
mean ocean density is the same as in the light-ocean model); vertical dotted lines mark the position 
of the mantle radius for the light-ocean and dense-ocean models. See Section 7.4 for details. 


thin dense layer), is constrained by the mantle density but not by the crust density. In 
the light-ocean model, 5p can thus be as large as 250kg/m^, which increases by 11%. 

The slope of the curve mainly depends on the relative radius of the mantle: it 
increases with xm so that it is slightly steeper for the dense-ocean model than for the 
light-ocean model. Fig. [T^ 


shows the change of (e^ given by Eq. (104)) as a function 
of the relative mantle radius for 6p = lOOkg/m^ (we known from Fig. 134 that /c2 depends 
linearly on dp so that it is easy to rescale the results for other 6p values). Curves have been 
drawn for = 0.59 (light-ocean model) but are nearly the same for = 0.67 (dense-ocean 
model). In the thin layer case (a = 0.01, dashed curve), the change of varies between 
0% and 15%. If a = 1, the change of k^ can be positive or negative depending on the 
balance between the mass concentration and screening effects. It is negative if xm ^ 3/4 
because the mass concentration effect wins over the screening effect. 

In conclusion, a thin and dense liquid layer at the bottom of a light ocean can increase 
the gravitational Love number by a significant amount (more than 10% for Titan) depend¬ 
ing on the radius of the mantle and on the density contrast with the core-mantle system. 
This is good news for the light-ocean models of Titan which otherwise predict that k 2 


should be much lower than what is measured with the Cassini spacecraft less et al. ,2012 
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8 Dynamical Love numbers 

8.1 Dynamical tides 

Tidal deformations of solid planetary bodies are most often computed in the static limit 
(equilibrium tides). Tides are indeed slow in comparison with seismic velocities so that 


inertial terms are negligible in the viscoelastic-gravitational equations (Section 2.3). By 


contrast, the fluid nature of stars and giant gaseous planets requires the inclusion of various 


Ogilvie 

2014 . The classical theory 

[e.g. Hendershott 

1981 . Dynamical 


of ocean tides on Earth is also intrinsically dynamical 
effects are thus potentially important in solid bodies with liquid layers, such as icy satellites 
with subsurface oceans. 

Dynamical tides, however, is a vast subject extending much beyond the viscoelastic- 
gravitational equations considered in this paper. First, rotation of the body leads to 
Coriolis and centrifugal effects which are of the same order of magnitude as inertial terms. 
Second, viscoelastic-gravitational equations do not incorporate a correct description of 
fluid dynamics (Navier-Stokes equations). These topics are clearly out of the scope of 
this paper. The Coriolis effect, for example, breaks spherical symmetry so that tidal 


deformations cannot be completely described with three scalar Love numbers (see Wahr 
1981 for the case of an oceanless Earth). I thus restrict myself to computing dynamical 


effects on the Love numbers of a spherically symmetric and non-rotating body enclosing 
a global fluid layer of zero viscosity. In this way, I hope to quantify the threshold beyond 
which dynamical effects become important for icy satellites. 

Dynamical corrections to the kn — hn formula are quantified by the parameter 
(Eq. (|74|) which depends on y^{Ri;)/yi{Ri;), the ratio of lateral to radial displacement 
of the fluid at the crust-ocean boundary. Computing requires solving the dynami¬ 
cal viscoelastic-gravitational problem within the ocean, but this cannot be done without 
solving the problem for the whole body at once. In the dynamical compressible case, 
the solution is obtained by numerically integrating the viscoelastic-gravitational equations. 
Numerical codes doing this job, however, usually diverge at low frequencies not far from 
the tidal frequencies of synchronously rotating satellites (note that the code included into 
SatStress is static). I will solve here the dynamical viscoelastic-gravitational problem for 


the rigid mantle model of Section 4.5 As before, inertial terms are neglected within solid 
layers. 

8.2 Dynamical incompressible liquid layer 


Denis et al. [T^ found a simple solution for dynamical tides in a homogeneous and 
incompressible fluid layer: 


= ar"-i+5r-’"-2. 
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2/3 

2/5 


j,n—l rf.—n—2 

a - b —— , 

n n + 1 


(111) 


The constants (a, b, a, (5) are fixed by the boundary conditions of the fluid layer. As seen in 
Section |2.41 vanishes in the fluid layer, while y 2 is linearly related to the other variables 


by Eq. (24). Finally, yg can be computed with Eq. The solution can be cast into the 


propagator matrix formalism: 


(2/1, 2/2 , 2/5 , yef = ^liq -{a, a,b, PY 


( 112 ) 


where the propagator matrix is given by Eq. (F.2) and its inverse by Eq. ( F.5[ ). Note 
that ys, instead of y2, is considered as a dependent variable because ya is discontinuous 


s. Using Eq. ( 

112 

) in the ocean 

Sabadini and 

Vermeersen 

2004 


one can solve for the tidal Love 


in the solid layers [e.g. 
numbers of an incompressible body made of an arbitrary number of homogeneous spherical 
layers. The ocean itself can be discretized in concentric layers of different densities. 


8.3 Rigid mantle model in its dynamical version 


The problem is further simplified by assuming an infinitely rigid mantle. In that case, the 
radial displacement of the mantle-ocean boundary vanishes: yi{Rm) = 0, where Rm is the 
mantle radius. This constraint reduces the number of unknown constants in Eq. (111). 
First, it yields a relation between a and b: 


b = —a R. 


2n+l 


(113) 


Combining Eqs. (Ill) and (113), I compute the ratio of tangential to radial displacements 
within the ocean: 

ya 1 n + 1 + n {Rm/rY'^'^^ 


2/1 


n (n-|-1) l — {Rm/rY^+^ 


(114) 


The assumption of an infinitely rigid mantle also imposes that /? = 0 in Eq. (Ill) 


This assertion can be proven by imposing the continuity of the variable yg (defined by 
Eq. (|17|)) at the mantle-ocean boundary. Intuitively, it can be understood as follows. In 
Eq. (Illll), the term of yg in is an interior gravity solution caused by the deformation of 


n 1 is 


the crust-ocean boundary and the layers above it. By contrast, the term of yg in r 
an exterior gravity solution caused by the deformation of the crust-mantle boundary and 
the layers beneath it. Thus the term of yg in vanishes if the mantle is infinitely rigid. 

Therefore y^ can be expressed within the ocean in terms of yg: 


/ n 

2/5 = - 2/5 

r 


(115) 


This relation is at the basis of the kn — proportionality derived in Appendix 
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8.4 Tilt factor 


The rigid mantle model, in its dynamical version, can be completely solved for the tilt 
factor and for the Love numbers. Substituting Eq. (114) into Eq. (74) yields a simple 
formula for the dynamical correction to the tilt factor: 


= -- 




n + 1 + 


(116) 


n(n + 1) 1 — 

where z = Rm/Re- The dynamical parameter is given in Table for Europa and Titan. 
If the ocean is shallow, 

-(‘ 1 ^) 

n(n + 1) D 

where D = Rg, — Rm is the ocean thickness. The main contributor to the tilt factor is 
typically the membrane spring constant Kt (Fig. 0. If the ocean is shallow, the 


relative magnitude of degree-two dynamical corrections in the tilt factor (Eqs. (@-(@) 

(118) 


IS 


Aoj 1 qio R 
X "" ~i2 "a m ■ 

For Europa, (/oj ~ 5 x 10“^ and fi = ft/(pgR) ~ 1.71 (if the crust is elastic, see Section 2.3) so 
that A[j/A ~ —2.4 x 10“®X/ {Dd). Fig. 14 shows how the tilt factor decreases as the ocean 
becomes shallower. The crust is elastic and either thick [d = 50 km) or thin [d = 10 km). 
If the crust is viscoelastic, the crust thickness should be interpreted as the lithospheric 
thickness: d —)■ dutho = \p/pE\d (Eq. (49)). The radius of the mantle varies with crust 
thickness and ocean depth. Otherwise, Europa is modeled as in Table (light density 
model) and Table If the ocean is deep (say D = 100 km), the dynamical correction 
decreases the tilt factor by a few percents: 1% if d = 50 km and 6% if d = 10 km. The 
dynamical correction becomes important if the ocean is less than 20 km thick; neglecting 
it leads to underestimate the crust thickness. The tilt factor can even become negative 
if the ocean is shallower than 1.2km (if d = 50km) or 5.9km (if d = 10km). All in all, 
dynamical corrections are significant for Europa unless the ocean is very deep. For Titan, 
dynamical corrections are much smaller because Qui is smaller by one order of magnitude. 

The rigid mantle model is based on rather restrictive assumptions. The hypothesis 
of an infinitely rigid mantle remains a good approximation for dynamical tides but can 
be dispensed with: the dynamical propagator matrix method is ideal for models with an 
arbitrary number of incompressible viscoelastic and liquid layers (Eq. 112). Ocean com¬ 


pressibility is expected to have a significant impact if the ocean is deep, but compressibility 
can be safely ignored in shallow oceans for which dynamical effects are maximum. 

8.5 Love numbers 

Without the static assumption, I cannot use gravity scaling in order to obtain explicit 
formulas for Love numbers. Instead, the kn — hn proportionality derived in Appendix 
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Figure 14: Dynamical corrections to the tilt factor of Europa as a function of ocean thickness. 
The tilt factor is divided by the relative crust thickness e = d/R. The asymptotic value is the tilt 
factor without dynamical corrections: 72 /e ~ 4.3 (similar to the elastic-like case of Fig. [^except 
that the imaginary part vanishes). See Section 8.4 for details. 


provides a supplementary constraint. Combining the kn — relation (Eq. with 


Eq. (G. 6 ), I get 


knr + 1 
hnr 


^nr + 1 


1 + 


1 


3g° _ 

272-|-1 

I tff, ^ 


ii^nr + 1 ) + Kpr) 


1 -|- (A -|- Atj) /l°^ -|- Ay. + Hpr 


(119) 


where (A:°^ -|- 1, /i°^) are the fluid-crust Love numbers if the fluid-crust density p° is equal 
to the ocean density p (implying 5p° = 5p). Terms beyond 0{e) are neglected. Density 
corrections are defined by 


pr - ^ ({n - 1 ) + 1 ) - ( 2 n + lfj — e, 


Kpr = 2 

Hpr = Kpr ~\~ Ap , 


( 120 ) 


where A' is defined by Eq. (G.5). Using the identity (C.2), one can check that {Kpr,H^ 


‘.prj 


are the rigid mantle limits of the density corrections {Kp, Hp) for the static case (Eq. (99)). 

If A(j = 0, the formulas for dynamical Love numbers (Eq. ( |119 )) are equivalent to the 
formulas for static Love numbers in the rigid mantle limit (Eq. (|98|)), the only difference 
being of Conversely, static formulas (for the rigid mantle model) are transformed 

into dynamical formulas by substituting 

A^A + A<^. (121) 

Being of different sign, the membrane spring constant A and the dynamical correction 


A^j have opposite effects (Fig. 14). The dynamical correction partially or totally cancels 
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the resistance of the crust to deformations, thus increasing Love numbers. What happens if 
the ocean thickness decreases even further? This question can be discussed with a simpler 
model in which the crust is incompressible and has the same density as the ocean. With 
these assumptions, it is unnecessary to require that ~ as done in Section 4.1 If 
= 0 and 6p = 0, the membrane formulas become 


{hnri hnr^ — 


1 


1 + (A + Atj) h° 


iKr,Kr)- 


( 122 ) 


This formula coincides with the thin shell limit of the dynamical homogeneous crust model 


with rigid mantle (Eqs. (E.6)-(E.7)), providing another check on the dynamical membrane 
formula. Eq. (122) shows that Love numbers diverge if A+A^^ = —l//i° 


resonance occurs. 


From Eq. (116), one sees that this cannot happen unless the ratio z = Rm/Re is close to 


one, i.e. the ocean must be shallow. The resonant ocean thickness is approximately equal 
to 

. (l -^e + ^e(A)] 

n(n + 1) V 2n + 1 ^ J 


D = 


-1 


(123) 


which is obtained by substituting Eqs. (117) and (C.l) into Eq. (122). I neglected terms of 


0(e) except the membrane spring constant which can be of 0(1) if the crust is hard (small 
or medium-sized satellite); other finite thickness corrections can be taken into account with 


the homogeneous crust model (Eq. (E.7)) but do not play a significant role here. 

In a large icy satellite, the crust is typically soft (A ~ 0(e), see Section 4.1) so that 


the resonant frequency is nearly the same as the one obtained by neglecting the crust 

(A^ = —1/L°^): 


= n (n + 1) 


1 — z 


2n+l 


+ 1 + 


n 


1 - 


2n + l^) R' 


(124) 


where z = Rm/R- This classical formula gives the frequency of the degree-n free oscillation 
of an incompressible fluid enclosing a rigid core (or mantle) \Lamb 1945, p. 454]. In astron¬ 
omy, this is called the / mode (/ for fundamental) or the surface gravity mode of oscillating 
stars and giant planets Ogilvie , 2014|. In seismology, one talks about the tsunami mode for 


obvious reasons Dahlen and Tromp ^ 1999 . In the model above (Eq. (122)), the resonance 


is damped by the viscoelasticity of the crust (proportional to Im(A)) but the damping is 
weak if the crust is soft. In more realistic models, other sources of damping arise such as 
solid/liquid friction or the viscosity of the fluid itself. Ogilvie [2014 illustrates the latter 
kind with the gravitational Love number of a homogeneous, incompressible, non-rotating, 
viscous fluid body: as in many classical oscillation problems, it amounts to the substitution 
—)■ -|- iaoj, where the coefficient a depends on the viscosity model (see his Eq. (16)). 

Fig. shows the enhancement of Europa’s Love numbers due to the dynamical res¬ 
onance. The crust is incompressible, is 10 km thick, and has the same density as the 
ocean (1000kg/m^). Otherwise, Europa is modeled as in Tables and(critical model). 
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Figure 15: Dynamical resonance enhancing the Love numbers of Europa for a shallow ocean. The 
solid, dashed, and dotted curves show the absolute value of / 12 , the absolute value of k 2 , and the 
imaginary part of ^ 2 , respectively. See Section 8.5|for details. 


Resonance occurs if the ocean is extremely shallow (D ~ 160 m), which is an unlikely 
configuration except in a transiting state: the ocean is either just born or on the point 
of freezing. At resonance, surface deformation (measured by I/ 12 I) and tidal dissipation 
within the crust (measured by Im{k 2 )) both diverge. If there is no dissipation in the core 


and mantle, tidal dissipation in the crust is related to surface deformation by Eq. (103): 
Im{k 2 ) = —(3/5)^|/i2p/m(A). Extremely large k 2 and /12 values should not be taken at 
face value but are rather a sign that the theory breaks down. In any case, viscoelastic- 
gravitational equations are only valid for small perturbations. Eurthermore, the present 
model does not take into account dissipation within the ocean and at its boundaries. 

Contrary to other dynamical approaches, the model presented here includes a crust 
above the ocean. However, it should not be considered as realistic for the analysis of 
dynamical resonances. Besides the lack of dissipation within the ocean, its fundamen¬ 
tal weakness is that other resonances appear once rotation is included, as is well-known 


from the analysis of Laplace tidal equations for a surface ocean [e.g. Chen et al ., 2014 


Matsuyama, 2014 Tyler, 2014 


9 Summary 

In a previous paper, I developed the massless membrane approach to tidal Love numbers of 
‘membrane worlds’, or icy satellites with subsurface oceans Beuthe. 2015 . This method. 


based on the classical equations of thin shell theory, leads to accurate formulas for Love 
numbers if the density contrast between the crust and the ocean is negligible. The error 
due to this approximation is never large because the ice and ocean densities do not differ by 
more than thirty percents, whereas viscoelasticity easily reduces the rigidity of the whole 
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crust by one order of magnitude. The membrane approach has the advantage of easily tak¬ 
ing into account crustal rheology through effective viscoelastic parameters. Nevertheless, 
the neglect of the crnst-ocean density contrast is an annoying featnre of this approach, 
especially now that measurements of Titan’s k 2 suggest that its ocean is very dense. Be¬ 
sides, it has a significant effect on the tilt factor which is useful quantity for crnst thickness 
estimates. Massless membrane theory is also unsatisfactory regarding the inclusion of crust 
compressibility because of its assumption that the npper and lower shell surfaces deform 
in the same way (Appendix [B|) . But let’s not throw ont the baby with the bath water: 
massless membrane theory remains largely valid thanks to its flexible formulation in terms 
of Love numbers. More precisely, the computation of Love numbers needs to be improved 
bnt the membrane formulas for tectonics and tidal dissipation are unchanged. 

For these reasons, I have rebuilt from scratch a membrane theory based on the viscoel¬ 
astic-gravitational theory used, among other things, for compnting tidal Love numbers in 
thick shell theory. The massive membrane approach starts with a perturbative expansion 
of the viscoelastic-gravitational equations in the small parameter e = d/R {d being the 
crust thickness): the variables of the problem are linearly propagated from the surface to 
the crust-ocean boundary. After integration over crust thickness, the viscoelastic response 
of the crust depends on three effective viscoelastic parameters fl, i>, and x (see Table for 
a list of parameters relevant to the membrane). These three complex numbers are directly 
computable for any depth-dependent rheology, thongh I benchmarked the formulas for a 
two-layer crust with Maxwell rheology. The compressibility factor x is new with respect to 
massless membrane theory and accounts for the compressibility problem mentioned above. 

At the crust-ocean boundary, two constraints (free slip and fluid constraint) lead to 
two relations between tidal Love numbers, here given at degree two: 


h 

k2 +I 


C- I - h2 , 

b + u 


6p 


(1 -I- At) /i2 - 5 — e, 
P 


(125) 


where dp is the crnst-ocean density contrast while Ay is the sum of the membrane spring 
constant A, the compressibility correction A^, the minor density correction Ap, and the 
dynamical correction A^^ (Eqs. ([70l)-([74l)). The I 2 —/i -2 and k 2 — ^2 relations are of 0{1) and 
0(e), respectively, and could in principle be computed to the next order of perturbation. 
These relations depend on the presence of a lithosphere (Eq. ©)■ There is thus no 
smooth transition, unless Sp = 0, between the kn — hr, relation and the hydrostatic relation 
fc°+l = /i° (Eq. @). 

In the static limit of equilibrium tides, one can write explicit formulas for Love numbers 
in terms of the fluid-crust Love numbers (fe^j ^ 2 )- fluid-crust model has the same internal 
strncture as the physical model except that the crust is fluid and of arbitrary density p° 
(in practice, p° is either equal to the ocean density or to the original crust density, see 


Section 5.2). The gravitational and radial Love numbers are given by Eq. (98) or, at 
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degree two, by 


^2 + 1 = 


h2 = 


fco + 1 


^ + 5 ( (^2 + 1) T+A + 


^ + (^ + i ^2) ^ 


(126) 


where ^° = p/pi is the ocean-to-bulk density ratio of the fluid-crust model (Eq. 

Mn ^ 


and 

{Kp, Hp) are density corrections proportional to 5p (Eq. (99)). These formulas are valid at 


0{e) if the crust is soft, as is generally the case for large icy satellites (Section |4.1[ ). If the 
crust is hard as could be the case for medium-sized and small icy satellites, A is of 0(1) 
instead of 0(e) and the formulas are only valid at 0 ( 1 ), meaning that density corrections 
should be neglected. 

Eluid-crust Love numbers are most easily computed with the propagator matrix ap¬ 
proach. The important thing to keep in mind is that the fluid-crust model does not have in 


general the same bulk density as the physical model (Eq. (85)). A useful two-layer model 
consists of a viscoelastic core-mantle below a homogeneous ocean (Eq. ( |C.3 )). Another sim¬ 
ple model is the three-layer body with an infinitely rigid mantle below two homogeneous 


fluid layers (Eq. (C.9)) which is a good toy model for ocean stratification. Nothing forbids. 


however, to compute fluid-crust Love numbers with more complicated models allowing for 
compressibility and continuous density variation. 

Benchmarking the membrane formulas with a numerical code shows that the predictions 
of k 2 , /i 2 and 72 have an error less than 1%, 0.5%, and 10%, respectively, assuming that 
the crust thickness is less than 5% of the surface radius. These error bounds correspond to 
worst-case scenarios in which the crust-ocean density contrast is maximum or the rheology 
of the crust is mostly fluid-like (Section]^. 

Besides their essential role in tidal deformations, tidal Love numbers serve to parameter¬ 
ize tidal heating. The membrane formulas for tidal Love numbers satisfy the micro-macro 
equivalence: the global heat flow from the whole body (proportional to Im{k 2 )) is the 
sum of the global heat flows from the crust (proportional to |/i 2 p/m(A)) and from the 
core-mantle (proportional to Im{k 2 )) (Section |5.6[ ). 

In a sense, the membrane approach factorizes the contributions from the crust and deep 
interior. The influence of the latter can thus be studied through fluid-crust Love numbers 
instead of physical Love numbers. In this way, one sees that the elasticity of the mantle 
increases the surface deformation by less than one percent, while the presence of a liquid 
core in Europa has an effect of less than two percents. It is thus a good approximation to 
consider the core-mantle system as being infinitely rigid (Sections |7.2[ ). This approximation 
is even better for smaller icy satellites because the nondimensional shear modulus of the 
mantle is inversely proportional to the surface radius and surface gravity (Eq. ( C.4| )). The 
rigid mantle approximation, however, is not essential. Eor example, the viscoelastic effect 
of a large liquid core can be taken into account in the two-layer fluid-crust model through 
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the equivalent shear modulus of the core-mantle system (Section |7.3[ ). In this paper, I 
emphasized the rigid mantle approximation in order to show that tidal effects in membrane 
worlds mainly depend on the shallow interior (crust and ocean). As long as crust thickness 
and rheology are not better known, it is besides the point to use sophisticated models of 
the deep interior to study tectonics and tidal dissipation. 

As is well-known, the deep interior structure affects Love numbers chiefly through 
global density stratification. For a hydrostatic body, Love numbers decrease if mass is 
more concentrated (as measured by the moment of inertia) but density stratification in 
an ocean sandwiched between elastic layers has a more complex effect. In particular, 
ocean stratification increases Love numbers if the mantle radius is large, especially if the 
ocean density abruptly increases near the bottom by large amount. This sereening effect 
occurs because the higher density layer at the bottom of the ocean screens the gravitational 
braking of the mantle. For Titan, a thin and dense liquid layer at the bottom of a light 
ocean can increase /c 2 by more than ten percents (Section |7.4[ ). 

The membrane approach is not limited to equilibrium tides. If the mantle is infinitely 
rigid and the ocean is homogeneous and incompressible, membrane formulas for dynamical 
Love numbers reveal the existence of a dynamical resonance for very shallow oceans (Sec¬ 
tion]^. Alternatively, dynamical Love numbers can be computed either with the dynamical 
homogeneous crust model (Appendix [E]), or more generally with the semi-dynamical prop¬ 
agator matrix method, where ‘semi’ means that dynamical effects are taken into account in 
the ocean but not in the solid layers (Appendix [F|) . On Europa, the resonance significantly 
decreases the tilt factor which becomes negative if the ocean is less than a few km thick (the 
exact threshold depends on the crust thickness). Dynamical effects should thus be taken 
into account when using the tilt factor for crust thickness estimates. More speculatively, 
the dynamical resonance could strongly increase tidal deformations and tidal heating in the 
crust if the ocean thickness is of the order of a few hundred meters, though the effect also 
depends on the existence of other resonances appearing once the rotation of the satellite is 
taken into account. 

Beyond tidal Love numbers, the membrane approach is well suited to the computation 
of load Love numbers which characterize the response of a spherically symmetric body 
to surface loading (Eq. Load Love numbers are related to tidal Love numbers by 


Eqs. (80)-(81). If the crust is incompressible and has a negligible density contrast with the 
ocean, the static load Love numbers of degree two are given by 


^2 

L' 


A /l2 — 1, 
AL 2 -5/(30 
1 + A 


(127) 


If the crust is thin and soft, these load Love numbers are close to ~ — 1 and h '2 - 5 /( 30 . 

If the crust is thick, these formulas can be transformed into thick shell formulas by the 
correspondence A o zuf (Appendix]^. This procedure bridges the gap between membrane 
formulas and the better-known Kelvin-Love formulas which are valid for a homogeneous 
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body. Load Love numbers of low degree are useful for planetary reorientation Matsuyamm 


et al. 

2014|, atmospheric loading Tokano et al. 

2011 , and ocean loading \Sohl et al. 

1995 

Tokano et al. 

2014 . 


For an easy implementation of the membrane formulas, I provide upon request the 
Mathematica notebook Membrane Worlds, nb which illustrates the membrane formulas with 
examples taken from this paper. 
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Appendix A Maxwell rheology 


The material in this section summarizes Appendix C of Beuthe 2015 . Under the assump¬ 


tion of zero bulk dissipation {Im{K)), viscoelastic parameters can be related to elastic 
constants, viscosity r] and angular frequency uj for any linear rheology. Once the depen¬ 
dence of the shear modulus /r on (rj, oj) is known, the viscoelastic Poisson’s ratio i' can be 
computed with 

He (1 + i^e) — M (1 “ 


E = 


‘^He (1 + i^e) + (1 — 21 /^) 


(A.l) 


in which the subscript E stands for ‘elastic’. This relation is valid for any linear rheology 
with zero bulk dissipation. Note that elastic incompressibility {ee = 1/2) implies that 
E = 1/2. 

For Maxwell rheology with no bulk dissipation, the viscoelastic shear modulus and 
Poisson’s ratio read 


He 

l-i 6 ’ 

3 ee -i{l + ee) S 
3 - 2 i{l + EE)d ’ 

where the dimensionless number 5 is defined by 

UJT] 


H = 

E = 


(A.2) 


(A.3) 


If the viscosity is high, ~ 0 so that h ^ He and e ^ ee- this is the elastic-like regime. If 
viscosity is low, 5 becomes large so that h ^ ^He/^ ^ and e ^ 1/2: this is the fluid-like 
regime (Fig. 1 of Beuthe 2015| ). 

The transition between the elastic- and fluid-like regimes occurs at <5 ~ 1, defining the 
critical state for which Im{g.) is maximum. Since the dissipated power per unit volume is 
proportional to Im{H), dissipation is nearly maximum in the critical state (‘nearly’, because 
the dissipated power also depends on /i 2 , see Eq. (101)). The critical viscosity associated 
with the critical state is given by 


Hcrit — 


HE 

U 


(A.4) 


Appendix B Massless membrane approach 


The massless membrane approach oi \Beuthe 2015] is based on the membrane theory of thin 
elastic shells. With respect to the massive membrane approach developed in the present 
paper, the massless membrane approach makes three supplementary approximations: 


1 . membrane of zero thickness. 
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2. no density contrast at the crust-ocean boundary, 

3. static approximation in all layers. 


The zero thickness assumption imposes that all terms of 0(e) are negligible in Eqs. (50)- 
(55) except those associated with the intrinsic properties of the membrane, which are its 
rigidity A and its surface mass density 6p x d (p is thus replaced by 6p). In this limit, 
Eqs. become (after substitution of Eqs. ([^-(f^) 


gyiiRe) 

-y2{Re) 

p 

yysiRs) 

y4{Rs) 

ybiRe) 

RyeiRe) 

RyviRe) 


hn, (B.l) 

— (A -|- Ap) hn e (2n 1) — , (B.2) 

P 

In , (B.3) 

0, (B.4) 

kn + l- 3e — hn, (B.5) 

Pb 


2n + I - 2,£— (r? - l) ^ hn , (B.6) 

Pb ' Xn + l + y 

2n + 1 — “ie — Khn — “ie — (2 (n — 1) hn — (2n + 1)). (B.7) 

Pb Pb 


One can show that these simpler membrane equations lead to the same formulas for Love 
numbers as the method assuming a massive membrane of finite thickness, except that the 
compressibility factor x is absent (A^^, = 0). For example, the kn — hn relation obtained 
with these equations is -|- 1 = (1 -|- A -|- Ap -|- A^)hn — (2n + l)(6p/p)e. Comparing 
Eq. (B.l) to Eq. (50), one sees that the absence of the compressibility term in thin shell 
theory results from the assumption that the upper and lower shell surfaces deform in the 
same way. 

If there is no density contrast at the crust-ocean boundary, the membrane must be 
massless {6p = 0), otherwise the membrane would carry a surface mass density in the limit 
of zero thickness. In this limit, Eqs. (B.1)-(B.6) tend to the membrane boundary conditions 
of the massless membrane approach: 


• Displacements are constant through the crust: yi(Re) = yi(R) and yz(Re) = y^(R)- 

• Gravity perturbations are constant through the crust: y^(Re) = yb{R) and y%(Re) = 

ye,(R)- 

• Stresses at the crust-ocean boundary are characterized by y^iRe) = 0 and y 2 (Re) = 
-pAgyi(R). 


To get a physical picture of the last equation, I rewrite it as 


q = (pyR) w , 


(B.8) 
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where U is the tidal potential, w = yi{Re)U is the membrane deflection and q = —y 2 {Re)U 
is the load acting on the bottom of the crust {Beuthe 2015 , Section 5.1). Eq. (B.8) has the 


form a Hooke’s law, justifying the name of membrane spring constant for A {Beuthe 2015 


Section 4.3). Eq. (B.2) shows how the Hooke’s law (Eq. ( B.8| )) of the massless membrane 
approach is modified by density corrections at 0{e). 

Einally the static approximation means that w = 0 and = 0. The three supplemen¬ 
tary assumptions of the massless membrane approach thus lead to A^ = Ap = A(^ = <5/3 = 0 
so that the kn — hn relation (Eq. (75)) rednces to kn + 1 = (1 + A)/i„, which is the kn — hn 
relation of the massless membrane approach Beuthe, 2015, Eq. (34)]. 


Appendix C Love numbers of fluid-crust models 


In Section |5.2[ the fluid-crust model is defined as having the same internal structure as 
the body with a viscoelastic crust (physical model), except that the crust is fluid and of 
density p°. The fluid-crust density p° is equal either to the original crust density {p° = p ), 
or to the density of the top layer of the ocean {p° = p). The latter choice makes the 
fluid-crust model simpler, but it changes the bulk density of the body from pi, to pi- I 
compute the fluid-crust Love numbers with the propagator matrix method [e.g. Sabadini 
and Vermeersen 2004 in the Fourier domain. This method is applicable if tides are static 


and the body is made of incompressible homogeneous layers. In the fluid layer, I use the 
variables ( 2 / 5 , y?) as defined by Saito [1974 . 


C.l. Rigid mantle and surface ocean 

If p° = p, the simplest fluid-crust model (next to the fluid body) is the two-layer incom¬ 
pressible body (bulk density pi) with an infinitely rigid mantle and a surface ocean (density 

-1 

, (C.l) 


/,o — 1 = 

'^nr ^nr ' 


2n + l 


where is the ocean-to-bulk density ratio for the fluid-crust model (Eq. (87)). This well- 
known formula is easily derived by combining the hydrostatic A:° — /i° relation (Eq. (lii)) 
with the proportionality relation (Eq. (G.l) in which ^ —)> ^°) as done for example in 


Section 4.5 of \Beuthe\ 2015| . Note that Eq. (C.l) does not depend explicitly on mantle 
parameters (radius and density): the mantle can be either large and light or small and 
dense. A useful identity is 

= + (C.2) 
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C.2. Viscoelastic mantle and surface ocean 


If p° = p, a more general fluid-crust model is the two-layer incompressible body (bulk 
density p^, radius R, surface gravity g) with viscoelastic core (radius Rm, shear modulus 
Pm) and surface ocean. The Love numbers of this model are given by 


^n{y,C) + Pm 

" B„(y,e) + (2n + l-3e)y^Arn’ 


where 






(C.3) 


(C.4) 


R pI Pb9°Ry 

where g° is the surface gravity of the fluid-crust model. If the mantle is infinitely rigid, 
Eq. (C.3) reduces to Eq. (C.l) as expected. The functions An and Bn are defined by 


^n{y,0 = fn{2n + l){l-^)pA{y,C,n), 

Bn{y,0 = fn{l- OPB{y,C,n) , (C.5) 


in which the superscript ° is omitted so as to simplify the notation. The factor /„ is given 
by 

fn = 


n 


2 (n - 1) (3 -h 4n 2n^) ’ 
while the functions pA and pB are defined by 

PA{y, C, n) = {2 {n-l) + 3y2^+^) (1 - 0 + (2n -M) ^, 

PB{y, C, n) = (2n + 1 - 30 (2 (n - 1) (1 - 0 + (2n + 1) y" c) - 9(1- 0 y'^+'C • 

The displacement of the mantle-ocean boundary is given by 


c = fn (2n + 1)' 


y 


,n+2 


( 1-0 


Bn{y, O) + (2n -I- 1 - 30) y^ Pm ’ 


(C.6) 


(C.7) 


(C.8) 


where /i™ = gyi{Rm)- For tides of degree two, these formulas are identical to Eqs. (118)- 
(120) of Beuthe\ [2015| . The various limits of this model (fluid or rigid core, uniform density, 
shallow ocean) are discussed in Appendix B of Beuthe [2015] . For example, /12 = 5/2 in the 
uniform density limit, in which case the mantle does not influence the surface deformation. 


C.3. Rigid mantle, ocean and fluid crust 


If p° = p, a simple fluid-crust model is the three-layer incompressible body (bulk density 
pD made of an infinitely rigid mantle, a homogeneous ocean (radius R — d, density p), and 
a homogeneous fluid crust (radius R, density p). The Love numbers of this model are 


+ 1 — (2^ + 1) 


yc(x,e°,r,n)-3(r-r) 




(2n + l- 3^°) pc{x, + 9C° (C - C) ’ 


(C.9) 


61 
















where 


and 


(x,r,f) 


R — d p p \ 


Pcix, C,C,n) = (2n + 1) (l - (l - x^) T ) + 3 (^ - T) . 


(C.IO) 

(C.ll) 


If the densities of the fluid crust and ocean are equal (^° = ^°) or if there is no fluid crust 
{x = 1), Eq. (C.9) reduces to Eq. (C.l) as expected. Similarly to Eq. (C.l), Eq. (C.9) does 
not depend explicitly on mantle parameters (radius and density). 

If the densities of the mantle and ocean are equal, the total mass equation yields 
^°x^ + ^°(1 — x^) = 1. Using this constraint to eliminate in favour of one can show 
that Eq. (C.9) reduces to An{x,^°)/Bn{x,^°), i.e. the formula for a two-layer fluid body 
(Eq. (C.3) with fim = 0). This result is expected because Love numbers do not depend on 
the properties of the mantle if there is no density contrast at the mantle-ocean boundary: 
this is the screening effect discussed in Section |7.4[ 


Appendix D Comparison with Wahr et al. 2006 


In Section |4.5[ I showed that the degree-two membrane formula for the tilt factor agrees 


with the analytical model of Wahr et al. 2006 in the rigid mantle limit. I will now do 


a similar comparison for the Love numbers themselves. Recall that the membrane Love 
numbers are expressed in terms of the fluid-crust Love numbers. If the fluid-crust density is 
equal to the ocean density, the fluid-crust model is an incompressible body of bulk density 
pI 7 ^ pfe made of an infinitely rigid mantle below a homogeneous ocean of density p reaching 
the surface (Eq. (C.l)). The fluid-crust Love numbers represent the first term, of 0{e^), 


in a perturbative expansion in the relative crust thickness s. 

also expand their solution about a two-layer model with infinitely 


Wahr et al. 2006 


rigid mantle and surface ocean. However, they choose a slightly different configuration for 
the density distribution at 0{£^). More precisely, the bulk density of their zeroth-order 
configuration is equal to the bulk density of the physical model (i.e. pb) instead of p^. This 
choice means the mantle density is implicitly adjusted in their zeroth-order configuration 
so as to compensate the density change due to the replacement of the crust by an ocean 


layer. The Love numbers of their zeroth-order configuration are thus given by Eq. (C.l) in 
which = p°/Pb is replaced by ^ = p/pb'. 


^0 + 1 — ^0 — ( 1 ~ 


3^ 


2n -|- 1 


-1 


(D.l) 


Using Eq. (85), I can express at 0(e) the Love numbers of the fluid-crust model in terms 
of Lq: 


fcnr + 1 = = ^0 ( 1 - (^ ( h) ^ £ 


-1 


(D.2) 
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Next, I substitute Eq. (D.2) into Eq. (98) and expand to 0{e) in the denominator. Reex¬ 


pressing the density terms with the help of the identity (C.2), I can write the Love numbers 
as 

/co + 1 


^nr “ 1 “ 1 — 


hnr — 


1 + 2I4T (t^ (^0 + 1) + Ro 
ho 


1 + A 1 + 


^0 ) + Ay -|- Hq 


where Kq and Hq are given by 

Ko 


Eor tides of degree two, 


Ho = 

Ko = 
Ho = 


3$ 

2n+l 


— {3ho + 2n 1) — £ , 

P 

Ap — (2n -2-1- 3ho) — £ ■ 
P 


.8-3C6p 

5-3^ p 

^ 40-3(3-z^)^ 5p ^ 


(D.3) 


(D.4) 


(5 z/) (5 - 30 P 

These equations generalize the results oi\Wahr et al. [2006] to a compressible crust with 


depth-dependent rheology. If the crust is incompressible {v = 1/2) and homogeneous 


{p = p), the expansions to 0(e) of Eq. (D.3) with {Ko, Hq) given by Eq. (D.5) agree with 
Eqs. (3)-(6) of Wahr et al. [2006[ . 


Appendix E Dynamical homogeneous crust model 

The homogeneous crust model describes a body in which (1) the crust is homogeneous and 
incompressible, and (2) there is a subsurface ocean, the top layer of which has the same 
density as the crust. Otherwise, the structure below the crust can be freely chosen. The 
static version of the model (for tides of degree two) is discussed in detail in Appendix A of 
Beuthe [2015 . I extend it here to dynamical tides. In comparison with static tides, there 
is mnch less freedom in choosing the deep interior structure because one cannot express 
the Love nnmbers in terms of arbitrary fluid-crust models. More complex models can be 
bnilt with the dynamical propagator matrix method (Appendix E). 

As elsewhere, {p, pb, g, R,d, p) denote the ocean (and crust) density, bulk density, sur¬ 
face gravity, surface radius, crust thickness and shear modulus of the crust, respectively. 
The Love numbers are expressed in terms of dimensionless ratios, three of which are 
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the fourth one being the dynamical correction (Eq. 

I obtain relations between Love numbers with the propagator matrix method as ex¬ 
plained in Appendix A of Beuthe [2015 : I propagate the y* from the surface down to the 


crust-ocean boundary, where I apply the free-slip condition (Eq. (23)) and the dynamical 


fluid constraint (Eq. (24)). The k 2 — /12 and I 2 — h 2 relations read 


k2 + l = {1 + Zhfj.)h2, 

h = zih2, 


(E.2) 


where 


Zh = 


zi = 


A^a jj- T Afe A.1^ 

p, -\- Ag A^ 
Ag p -|- Xd 

p T Ag A^ 


(E.3) 


in which Xj are geometrical factors depending only on x (Table [TT|. 
Load Love numbers satisfy similar relations: 


k'2 + 1 — (1 -|-Z/j/i)/12-b — , 

^2 ~ Zih 2 , 


(E.4) 


The gravitational Load number is related to the tidal Love numbers by the Saito-Molodensky 


relation (Eq. (80)) which can be combined with the k 2 — /12 and k '2 — h '2 relations in order 
to express (/c^j h' 2 ) in terms of / 12 . Load and tidal Love numbers are thus related by 


k '2 = Zhph2 -I, 


h '2 = 


Zhph2- 5/(3g) 

1 + Zhp 


(E.5) 


Alternatively, these equations can be obtained from the membrane formulas (Eq. (127)) 
by the substitution A —)■ Zhp, as already noted in Beuthe 2015| (see his Table 8 and 
Appendix A). 

In the thin shell limit, the factors appearing in the Love number relations behave to 
leading order as 


Zhp ~ 

Zl ~ 


A -|- A(^ , 

11 V 6/i 


(E.6) 


where A 


(24/ll)/ie is the incompressible membrane spring constant (Eq. (0). 
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Table 11: Geometrical factors Xj appearing in the Love number relations (Eq. (E.31). In each case, 
the polynomial in the denominator is P = 24 + 40a;^ — 45 — 19x^°. The polynomials appearing 

below have a simple or double root at x = 1, but factorizing them yields longer expressions. ‘Thin 
shell’ denotes the thin shell limit (e = 1 — x = d/R)\ the coefficients of fi and are expanded to 
next-to-leading and leading orders in e, respectively. ‘Homog.’ denotes the homogeneous limit in 
which the shell hlls the whole body, i.e. x = 0. The factors Xa and X^. are shown as functions of x 
in Fig. 13 of . 


Beuthe 


2015 



Thick shell 

Thin shell 

Homog. 

Aa 

f (19 - 75 x^ -M12 x^ - 75 x^ -h 19 x^^) /P 


19 

5 

Ab 

(19 -h 45x3 - 40x^ - 24x3°) /P 

1 

0 

Ac 

1 (36 - 100x3 -t 308x° - 225 - 19x3°) /P 

11 V 33 / 

3 

10 

Ad 

1 X^ (3 -1- 5x3 _ _|_ 2a;10) /P 

1 

22 

0 

Ac 

5 X^ (1 — x3 — x"^ -1- x3°) /P 

11 ^ 

0 


Some assumptions about the deep interior are required in order to determine . In the 
rigid mantle model, is given by Eq. (116)andA:2 is proportional to ^2 (Eq. (G.l)). Com¬ 
bining the latter constraint with the /c 2 — /12 relation, I get explicit formulas for dynamical 
Love numbers: 


(^2 ) ^-2 , ^2 , ^ 2 ) — 


1 


1 + A 


h° h° —1 _ 

'^2r 1 '^2r 1 ^ 1 


(E.7) 


where + 1 = /i 2 ^ = 5/(5 — 3^) (Eq. (C.l)). Alternatively, one can derive this formula 
with the dynamical propagator matrix method (Appendix F). 

Additionally, if tides are static and the body is of uniform density, Zh = Xa and 
/i 2 = A :2 + 1 = 5/2. Such a model describes a two-layer body made of a liquid core 
surrounded by a solid layer: 


(/c2 ,h2 jk'o , ho) = -F- (- , - — 1, —^ . 

^ l + ^Xafi V 2 2 ’ sy 


(E.8) 


This equation reduces to the well-known Kelvin-Love formula if the solid layer extends to 
the center of the body {Xa = 19/5 if x = 0): 


(^2 ) ^2 , ^2 , h ^) — 


1 


l + fjl \2’2 


3 5 


- 1 ,- 


(E.9) 


as quoted by Lambeck |1980| (his Eq. (2.1.9)). 
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Appendix F Dynamical propagator matrix 


In the propagator matrix method, one usually assumes that tides are static and the body 
is incompressible. If either of these assumptions does not hold, the elastic-gravitational 
problem can be solved in terms of spherical Bessel functions if grjT is constant within 
the body, Qr being the gravitational acceleration at radius r. \Love ](mT 


and Pekeris 


and Jarosch 1958 computed in this way the deformations of a compressible body of 


uniform density. Assuming that gr/r is approximately constant inside each layer, Gilbert 


and Backus 1968| proposed a propagator matrix method, in the Fourier domain, for the 
oscillations and dynamical tides of an elastic stratified compressible body. | Vermeersen et al. 
later considered the viscoelastic problem in the Laplace domain in order to model 


1996 


deformations at the geological time scale (see also Appendix A of Sabadini and Vermeersen 


2004|). This method has not become popular for several reasons. First, the approximation 


of constant gr/r in each layer introduces errors that possibly offset any gain due to the 
inclusion of compressibility, at least in models with few layers (benchmarking has yet to 
be done, as observed by Vermeersen et al. 1996] ). Second, spherical Bessel functions are 
a bit tricky to use, the more so in the viscoelastic generalization of the problem. Third, 
Bessel functions obscure the dependence of Love numbers on physical parameters, so that 
the analytical formulas cannot be interpreted without being first numerically evaluated. 
Fourth, the static and incompressible limits are problematic so that comparisons with 
simpler cases are difficult. In conclusion, the benefits brought by this rather complicated 
analytical method are not obvious: you can as well resort to numerical integration. 

By contrast, the dynamical solution within a homogeneous incompressible liquid layer 
takes a very simple power form (Section |8.2| ). Since dynamical effects are much smaller in 
solid than in liquid layers, I propose to use the static and dynamical solutions within the 
solid and liquid layers, respectively. The static propagator matrix appears in the literature 
(e.g. Sabadini and Vermeersen 2004] with other m conventions) but the dynamical prop¬ 
agator matrix for an incompressible liquid layer does not seem to have been written down 
before. This propagator matrix does not involve the variables ys and y^: the former is dis¬ 
continuous at the fluid-solid interfaces and is considered as a dependent variable (through 
Eq. (24)) whereas the latter vanishes everywhere within the fluid. Combining Eq. ( jlll ) 
with Eqs. (18) and (24) gives the following propagation formula: 


(yi ,y2,y5, yaf = ^Uq ■ {a, a, b, , (E.l) 

where (a, a, 6, /3) are free constants. In the conventions of Sabadini and Vermeersen [2004] , 
the left-hand side is replaced by (yi, 2/3, —ys, —ye)^- ^Hq is the 4x4 propagator matrix 
defined by 

Yliq = Y, 

where Cuq is a diagonal matrix given by 


Ciiq = diag (r 


n—1 



(E.2) 


(E.3) 
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while is given by 


( ^ 

0 

1 

0 

\ 

p (ft - 

-P 

p (ft +st) 

-P 


0 

1 

0 

1 


\ —AirGp 

2n+l 

r 

—AnGp 

0 

/ 


(F.4) 


where a; is the tidal angular frequency, p is the uniform density of the layer, and pr is the 
unperturbed gravitational acceleration at radius r. 

The inverse propagator matrix reads 


■VT—1 _ -y 

^ liq ~ ‘ ^ Itq ^ 


(F.5) 


where Duq is a diagonal matrix given by 


D 


liq 


2n + l 
while Y/jg is given by 


* diag r- . 


puj^ 


puj^ 


(F. 6 ) 


Yuq = 


^ P{sr + ^,} -1 0 

AttGp 

-p (ft - 

Y —AnGp 


\ 


0 0 1 

1 p 0 

0 ^ -1 / 


(F.7) 


As an illustration, consider a three-layer body with viscoelastic core, ocean and crust. The 
viscoelastic-gravitational problems depends on 13 free constants (3 in the core, 4 in the 
ocean and 6 in the crust). Besides the 3 surface boundary conditions (Eq. ®) , there are 8 
constraints due to the continuity of the variables (yi, 2 / 2 ) ?/ 5 , ye) at the core-ocean and crust- 
ocean boundaries and 2 constraints due to 2/4 = 0 at the same boundaries. The problem 
can thus be solved by propagating the general solution from the center to the surface and 
applying the various boundary conditions [e.g. Roberts and Nimmo. 2008 . If the ocean 


density increases with depth, one can discretize the ocean in layers of uniform density. 


Appendix G kn — hn proportionality 

The rigid mantle model is a three-layer body made of an infinitely rigid mantle, a homo¬ 
geneous and incompressible ocean, and a crust. In this model, the deep interior does not 
contribute to the gravity perturbation: kn is only due to the crust deformation and to the 
ocean bulge just below the crust. If the crust is incompressible and has the same density 
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p as the ocean, the induced geoid perturbation is proportional to the surface deformation 

3? 


Beuthe, 2015 


k — 

i^n.r — 


hr, 


2n + 1 

where ^ = p/pb as before and the subscript r denotes the rigid mantle model. 


(G.l) 


What is the analogue of Eq. (G.l) if the crust is compressible and does not have the 


same density as the ocean? Although there is no similar relation for a thick crust, I will 


show how to extend Eq. (G.l) if the crust is thin. For this purpose, I evaluate yg in terms 


of (yi, ys) on the ocean side of the crust-ocean boundary with Eq. © in which I substitute 
Eq. (flT^: 


^ ~ ^9yi{Re) • 


(G.2) 


Since (yi,y 5 ,y 6 ) are continuous across the crust-ocean boundary, I can replace these vari¬ 


ables by their values on the crust side of the boundary (Eqs. (50)-(55)), before substituting 


Eqs. (67)-(68) for Sn and Expanding to 0(e), I obtain the analog of Eq. (G.l) for a 


massive membrane of finite thickness: 


hnr — 


3? 


2n -|- 1 


1 Ay “t“ An “t“ 3 — E I h' 


(G.3) 


where Ay and Ap are defined by Eqs. (72)-(73) 


As a last step, I express the kn — hn proportionality in terms of the Love numbers of 


the fluid-crust model, with p° = p (implying 6p° = 6p). Using Eq. (85), I rewrite Eq. (G.3) 


so that it depends on = p/p^ instead of ^ = p/pb- This substitution is only necessary 
in the term of 0(1) because ^ and are interchangeable in terms of 0(e). The kn — hn 
proportionality becomes 


knr — 


3e 


2n -|- 


- (l -F Ay -F A'p) hn 


where 


R'p — Rp A ^ (1 — C) ^ E . 


Finally I substitute Eq. (G.l) into Eq. (G.4) so that h'. 


(G.4) 

(G.5) 

appears explicitly. In its final 


form, the kn — hn proportionality reads 


knr — (l — (hnr) (^ + ^r. 


(G.6) 


Though the induced geoid perturbation is still proportional to the radial displacement. 


Eqs. (G.l) and (G.6) differ by two types of terms: (1) Ay is the compressibility correction 
appearing when the membrane is of finite thickness, and (2) A^ is the correction due to 
the density contrast at the crust-ocean boundary. 
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